

### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA




prepanel.default.bwplot <-
    function(x, y, box.ratio,
             horizontal = TRUE,
             levels.fos = length(unique(y)), ...)
{
    temp <- .5  #* box.ratio/(box.ratio+1)
    if (horizontal)
        list(xlim = range(x[is.finite(x)]),
             ylim = c(1-temp, levels.fos + temp),
             dx = 1,
             dy = 1)
    else 
        list(xlim = c(1-temp, levels.fos + temp),
             ylim = range(y[is.finite(y)]),
             dx = 1,
             dy = 1)
}





panel.barchart <-
    function(x, y, box.ratio=1, horizontal = TRUE, col=bar.fill$col, ...)
{
    bar.fill <- trellis.par.get("bar.fill")
    if (horizontal) {
        xmin <- current.viewport()$xscale[1]
        height <- box.ratio/(1+box.ratio)
        
        for (i in seq(along=x)) {
            grid.rect(gp = gpar(fill = col),
                      y = y[i],
                      x = unit(0,"npc"),
                      height = height,
                      width = x[i]-xmin,
                      just = c("left", "centre"),
                      default.units = "native")
        }
    }
    else {
        ymin <- current.viewport()$yscale[1]
        width <- box.ratio/(1+box.ratio)
        
        for (i in seq(along=y)) {
            grid.rect(gp = gpar(fill = col),
                      x = x[i],
                      y = unit(0,"npc"),
                      height = y[i] - ymin,
                      width = width,
                      just = c("centre", "bottom"),
                      default.units = "native")
        }
    }
}






panel.dotplot <-
    function(x, y, horizontal = TRUE,
             pch = dot.symbol$pch,
             col = dot.symbol$col,
             lty = dot.line$lty,
             lwd = dot.line$lwd,
             col.line = dot.line$col,
             levels.fos = NULL, ...)
{
    dot.line <- trellis.par.get("dot.line")
    dot.symbol <- trellis.par.get("dot.symbol")
    if (horizontal) {
        yscale <- current.viewport()$yscale
        if (is.null(levels.fos)) levels.fos <- floor(yscale[2])-ceiling(yscale[1])+1
        panel.abline(h=1:levels.fos, col=col.line,
                     lty=lty, lwd=lwd)
        panel.xyplot(x = x, y = y, col = col, pch = pch, ...)
    }
    else {
        xscale <- current.viewport()$xscale
        if (is.null(levels.fos)) levels.fos <- floor(xscale[2])-ceiling(xscale[1])+1
        panel.abline(v=1:levels.fos, col=col.line,
                     lty=lty, lwd=lwd)
        panel.xyplot(x = x, y = y, col = col, pch = pch, ...)
    }
}





panel.stripplot <-
    function(x, y, jitter.data = FALSE, factor = 0.5, horizontal = TRUE,
             pch=plot.symbol$pch, col=plot.symbol$col, ...)
{

    plot.symbol <- trellis.par.get("plot.symbol")
    y.jitter  <-
        if (horizontal && jitter.data) jitter(y, factor = factor)
        else y
    x.jitter  <-
        if (!horizontal && jitter.data) jitter(x, factor = factor)
        else x
    panel.xyplot(x = x.jitter, y = y.jitter, pch = pch, col = col, ...)
    
}





panel.bwplot <-
    function(x, y, box.ratio=1, horizontal = TRUE, pch=box.dot$pch,
             col=box.dot$col, cex = box.dot$cex,
             levels.fos = NULL, ...)
{
    
    box.dot <- trellis.par.get("box.dot")
    box.rectangle <- trellis.par.get("box.rectangle")
    box.umbrella <- trellis.par.get("box.umbrella")
    plot.symbol <- trellis.par.get("plot.symbol")

    ## In case levels.fos is not given (which should not happen), I'll
    ## be working on the premise that EACH INTEGER in the y-RANGE is a
    ## potential location of a boxplot. 

    if (horizontal) {
        yscale <- current.viewport()$yscale
        if (is.null(levels.fos)) levels.fos <- floor(yscale[2])-ceiling(yscale[1])+1
        lower <- ceiling(yscale[1])
        height <- box.ratio/(1+box.ratio)
        xscale <- current.viewport()$xscale
        if (levels.fos > 0)
            for (i in 1:levels.fos) {

                yval  <- i
                stats <- boxplot.stats(x[y==yval])
                
                
                if (stats$n>0)
                {
                    push.viewport(viewport(y=unit(yval, "native"),
                                           height = unit(height, "native"),
                                           xscale = xscale))
                    
                    r.x <- (stats$stats[2]+stats$stats[4])/2
                    r.w <- stats$stats[4]-stats$stats[2]
                    grid.rect(x = unit(r.x, "native"), width = unit(r.w, "native"),
                              gp = gpar(lwd = box.rectangle$lwd,
                              lty = box.rectangle$lty,
                              col = box.rectangle$col))
                
                    grid.lines(x = unit(stats$stats[1:2],"native"),
                               y=unit(c(.5,.5), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.lines(x = unit(stats$stats[4:5],"native"),
                               y=unit(c(.5,.5), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.lines(x = unit(rep(stats$stats[1],2),"native"),
                               y=unit(c(0,1), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.lines(x = unit(rep(stats$stats[5],2),"native"),
                               y=unit(c(0,1), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.points(x=stats$stats[3], y=.5, pch=pch, 
                                size = unit(cex * 2.5, "mm"),
                                gp = gpar(col = col, cex = cex))
                    
                    if ((l<-length(stats$out))>0)
                        grid.points(x = stats$out, y = rep(.5,l),
                                    size = unit(plot.symbol$cex * 2.5, "mm"),
                                    pch = plot.symbol$pch,
                                    gp = gpar(col = plot.symbol$col,
                                    cex = plot.symbol$cex))
                    
                    pop.viewport()
                
                }
            }
    
    }
    else {
        xscale <- current.viewport()$xscale
        if (is.null(levels.fos)) levels.fos <- floor(xscale[2])-ceiling(xscale[1])+1
        lower <- ceiling(xscale[1])
        width <- box.ratio/(1+box.ratio)
        yscale <- current.viewport()$yscale
        if (levels.fos > 0)
            for (i in 1:levels.fos) {
                xval  <- i
                stats <- boxplot.stats(y[x==xval])

                if (stats$n>0)
                {
                    push.viewport(viewport(x = unit(xval, "native"),
                                           width = unit(width, "native"),
                                           yscale = yscale))
                    
                    r.x <- (stats$stats[2]+stats$stats[4])/2
                    r.w <- stats$stats[4]-stats$stats[2]
                    grid.rect(y = unit(r.x, "native"), height = unit(r.w, "native"),
                              gp = gpar(lwd = box.rectangle$lwd,
                              lty = box.rectangle$lty,
                              col = box.rectangle$col))
                
                    grid.lines(y = unit(stats$stats[1:2],"native"),
                               x = unit(c(.5,.5), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.lines(y = unit(stats$stats[4:5],"native"),
                               x = unit(c(.5,.5), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.lines(y = unit(rep(stats$stats[1],2),"native"),
                               x = unit(c(0,1), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.lines(y = unit(rep(stats$stats[5],2),"native"),
                               x = unit(c(0,1), "npc"),
                               gp = gpar(col = box.umbrella$col,
                               lwd = box.umbrella$lwd,
                               lty = box.umbrella$lty))
                    
                    grid.points(y = stats$stats[3], x = .5, pch = pch, 
                                size = unit(cex * 2.5, "mm"),
                                gp = gpar(col = col, cex = cex))
                    
                    if ((l<-length(stats$out))>0)
                        grid.points(y = stats$out, x = rep(.5,l),
                                    size = unit(plot.symbol$cex * 2.5, "mm"),
                                    pch = plot.symbol$pch,
                                    gp = gpar(col = plot.symbol$col,
                                    cex = plot.symbol$cex))
                    
                    pop.viewport()
                
                }
            }
    
    }
}



dotplot <-
    function(formula,
             data = parent.frame(),
             panel = "panel.dotplot",
             prepanel = NULL,
             strip = TRUE,
             groups = NULL,
             horizontal = NULL,
             ...,
             subset = TRUE)
{

    ## m <- match.call(expand.dots = FALSE)
    dots <- list(...)
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    formname <- deparse(substitute(formula))
    form <- eval(substitute(formula), data, parent.frame())

    call.list <- c(list(groups = groups, subset = subset, horizontal = horizontal,
                        panel = panel, prepanel = prepanel, strip = strip,
                        box.ratio = 0),
                   dots)

    if (!inherits(form, "formula") && is.numeric(form)) {
        call.list$formula <- form
        if (is.null(call.list$xlab)) call.list$xlab <- formname
        else if (is.list(call.list$xlab) && is.null(call.list$xlab$lab))
            call.list$xlab$label <- formname 
    }
    else {
        call.list$formula <- formula
        call.list$data <- data
    }

    do.call("bwplot", call.list)

    ## lapply(dots, eval, data, parent.frame())))
}



barchart <-
    function(formula,
             data = parent.frame(),
             panel = "panel.barchart",
             prepanel = NULL,
             strip = TRUE,
             box.ratio = 2,
             groups = NULL,
             horizontal = NULL,
             ...,
             subset = TRUE)
{

    ## m <- match.call(expand.dots = FALSE)
    dots <- list(...)
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    do.call("bwplot",
            c(list(formula = formula, data = data, horizontal = horizontal,
                   groups = groups, subset = subset,
                   panel = panel, prepanel = prepanel, strip = strip,
                   box.ratio = box.ratio),
              dots))
    ## lapply(dots, eval, data, parent.frame())))
}


stripplot <-
    function(formula,
             data = parent.frame(),
             panel = "panel.stripplot",
             prepanel = NULL,
             strip = TRUE,
             jitter = FALSE,
             factor = .5,
             box.ratio = if (jitter) 1 else 0,
             groups = NULL,
             horizontal = NULL,
             ...,
             subset = TRUE)
{
    ## m <- match.call(expand.dots = FALSE)
    dots <- list(...)
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    formname <- deparse(substitute(formula))
    form <- eval(substitute(formula), data, parent.frame())

    call.list <- c(list(groups = groups, subset = subset, horizontal = horizontal,
                        panel = panel, prepanel = prepanel, strip = strip,
                        jitter = jitter, factor = factor,
                        box.ratio = 0),
                   dots)
    
    if (!inherits(form, "formula") && is.numeric(form)) {
        call.list$formula <- form
        if (is.null(call.list$xlab)) call.list$xlab <- formname
        else if (is.list(call.list$xlab) && is.null(call.list$xlab$lab))
            call.list$xlab$label <- formname 
    }
    else {
        call.list$formula <- formula
        call.list$data <- data
    }

    do.call("bwplot", call.list)

    ## lapply(dots, eval, data, parent.frame())))
}


bwplot <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.bwplot",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             box.ratio = 1,
             horizontal = NULL,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ##m <- match.call(expand.dots = FALSE)
    ##dots <- m$...
    ##dots <- lapply(dots, eval, data, parent.frame())

    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing

    formname <- deparse(substitute(formula))
    formula <- eval(substitute(formula), data, parent.frame())

    form <-
        if (inherits(formula, "formula"))
            latticeParseFormula(formula, data)
        else {
            if (!is.numeric(formula)) stop("invalid formula")
            else {
                list(left = rep("", length(formula)),
                     right = formula,
                     condition = NULL,
                     left.name = "",
                     right.name = formname)
            }
        }
    if (is.null(form$left)) form$left <- rep("", length(form$right))

    cond <- form$condition
    number.of.cond <- length(cond)
    y <- form$left
    x <- form$right

    if (is.null(horizontal)) {
        horizontal <-
            if ((is.factor(x) || is.shingle(x)) && is.numeric(y)) FALSE
            else TRUE
    }

    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())

    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    if (subscripts) subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    y <- y[subset, drop = TRUE]
    if (subscripts) subscr <- subscr[subset, drop = TRUE]

    if (horizontal) {
        if (!(is.numeric(x))) {
            x <- as.numeric(x)
            warning("x should be numeric")
        }
        y <- as.factorOrShingle(y)
        is.f.y <- is.factor(y)  # used throughout the rest of the code
        num.l.y <- nlevels(y)

        if (missing(xlab)) xlab <- form$right.name
        if (missing(ylab)) ylab <- if (is.f.y) NULL else form$left.name
    }
    else {
        if (!(is.numeric(y))) {
            y <- as.numeric(y)
            warning("y should be numeric")
        }
        x <- as.factorOrShingle(x)
        is.f.x <- is.factor(x)  # used throughout the rest of the code
        num.l.x <- nlevels(x)

        if (missing(ylab)) ylab <- form$left.name
        if (missing(xlab)) xlab <- if (is.f.x) NULL else form$right.name
    }

    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- form$right.name
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <- form$left.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    ##scales <- 
    ##if (missing(scales)) scales 
    ##else eval(m$scales, data, parent.frame())


    ## The following is to make the default alternating FALSE for factors
    if (is.character(scales)) scales <- list(relation = scales)
    if (is.null(scales$alternating)) {
        if (horizontal) {
            if (is.null(scales$y)) scales$y <- list(alternating = FALSE)
            else if (is.null(scales$y$alternating)) scales$y$alternating <- FALSE
        ## bug if y="free" ? but who cares
        }
        else {
            if (is.null(scales$x)) scales$x <- list(alternating = FALSE)
            else if (is.null(scales$x$alternating)) scales$x$alternating <- FALSE
        ## bug if x="free" ? but who cares
        }
    }
    foo <- c(foo,
             do.call("construct.scales", scales))

    if (horizontal) {
        if (is.logical(foo$y.scales$at)) foo$y.scales$at <- 1:num.l.y
        if (is.f.y && is.logical(foo$y.scales$labels))
            foo$y.scales$labels <- levels(y)
    }
    else {
        if (is.logical(foo$x.scales$at)) foo$x.scales$at <- 1:num.l.x
        if (is.f.x && is.logical(foo$x.scales$labels))
            foo$x.scales$labels <- levels(x)
    }

    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        ## warning("Are you sure you want log scale for y ?")
        ylog <- foo$y.scales$log
        ybase <-
            if (is.logical(ylog)) 10
            else if (is.numeric(ylog)) ylog
            else if (ylog == "e") exp(1)

        y <- log(y, ybase)
        if (have.ylim) ylim <- log(ylim, ybase)
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)|is.na(y)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    foo$panel.args.common <- dots
    foo$panel.args.common$box.ratio <- box.ratio
    foo$panel.args.common$horizontal <- horizontal
    foo$panel.args.common$levels.fos <- ## fos - the factor/shingle in x/y
        if (horizontal) num.l.y else num.l.x
    if (subscripts) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    if (horizontal) {
                        if (is.f.y) 
                            foo$panel.args[[panel.number]] <-
                                list(x = x[id],
                                     y = as.numeric(y[id]))
                            
                        else {  # shingle
                            panel.x <- numeric(0)
                            panel.y <- numeric(0)
                            for (k in 1:num.l.y) {
                                tid <- id & (y >= levels(y)[[k]][1]) & (y <= levels(y)[[k]][2])
                                panel.x <- c(panel.x, x[tid])
                                panel.y <- c(panel.y, rep(k,length(tid[tid])))
                            }
                            foo$panel.args[[panel.number]] <-
                                list(x = panel.x,
                                     y = panel.y)
                        }
                    }
                    else {
                        if (is.f.x)
                            foo$panel.args[[panel.number]] <-
                                list(x = as.numeric(x[id]),
                                     y = y[id])
                            
                        else {  # shingle
                            panel.x <- numeric(0)
                            panel.y <- numeric(0)
                            for (k in 1:num.l.x) {
                                tid <- id & (x >= levels(x)[[k]][1]) & (x <= levels(x)[[k]][2])
                                panel.y <- c(panel.y, y[tid])
                                panel.x <- c(panel.x, rep(k,length(tid[tid])))
                            }
                            foo$panel.args[[panel.number]] <-
                                list(x = panel.x,
                                     y = panel.y)
                        }
                    }


                    if (subscripts)
                        foo$panel.args[[panel.number]]$subscripts <-
                            subscr[id]

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.bwplot,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}



















### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA




prepanel.default.cloud <-
    function(distance, xlim, ylim, zlim, zoom = 1,
             rot.mat = rot.mat, aspect = aspect, ...)
{
    aspect <- rep(aspect, length=2)

    corners <-
        rbind(x = c(-1,1,1,-1,-1,1,1,-1) / 2,
              y = c(-1,-1,-1,-1,1,1,1,1) / 2 * aspect[1],
              z = c(-1,-1,1,1,-1,-1,1,1) / 2 * aspect[2])

    ## box ranges and lengths
    cmin <- lapply(corners, min)
    cmax <- lapply(corners, max)
    clen <- lapply(corners, function(x) diff(range(x)))

    ##if (!missing(R.mat))
    ##  {;} #do something

    corners <- rot.mat %*% corners

    zback <- min(corners[3,])
    zfront <- max(corners[3,])
    za <- (zfront * (1-distance) - zback) / (zfront - zback)
    zb <- distance / (zfront - zback)

    corners[1,] <- (za + zb * corners[3,]) * corners[1,]
    corners[2,] <- (za + zb * corners[3,]) * corners[2,]

    xrng <- range(corners[1,])
    yrng <- range(corners[2,])
    slicelen <- max(diff(xrng), diff(yrng))

    list(xlim = extend.limits(xrng, length = slicelen) / zoom,
         ylim = extend.limits(yrng, length = slicelen) / zoom,
         dx = 1, dy = 1)
}











panel.cloud <-
    function(x, y, z, subscripts, distance, xlim, ylim, zlim,
             subpanel = "panel.xyplot",
             rot.mat = rot.mat, aspect = aspect,
             zcol, col.regions, par.box = NULL,
             ## next few arguments are an attempt to support
             ## scales. The main problem with scales is that it is
             ## difficult to figure out the best way to place the
             ## scales. Here, they would need to be specified
             ## explicitly. Maybe this code can be used later for a
             ## proper implementation
             xlab, ylab, zlab, scales.3d,
             proportion = 0.6, wireframe = FALSE,
             scpos = list(x = 1, y = 8, z = 12), 
             ...)
{
    if (any(subscripts)) {

        par.box.final <- trellis.par.get("box.3d")
        if (!is.null(par.box)) par.box.final[names(par.box)] <- par.box

        subpanel <-
            if (is.character(subpanel)) get(subpanel)
            else eval(subpanel)

        aspect <- rep(aspect, length=2)

        x <- x[subscripts]
        y <- y[subscripts]
        z <- z[subscripts]

        if (wireframe) {
            ord <- order(y, x)
            x <- x[ord]
            y <- y[ord]
            z <- z[ord]
            nx <- length(unique(x))
            ny <- length(unique(y))
            len <- length(z)
            if (nx * ny != len) stop("Incorrect arguments")
            zcol <- rep(zcol, len)

            idtf <- c(rep( c(rep(TRUE, nx-1), FALSE), ny-1 ), rep(FALSE, nx))
            id0 <- (1:(nx*ny))[idtf]
            id1 <- id0 + 1
            id2 <- id1 + nx
            id3 <- id0 + nx
        }


        corners <-
            data.frame(x = c(-1,1,1,-1,-1,1,1,-1) / 2,
                       y = c(-1,-1,-1,-1,1,1,1,1) / 2 * aspect[1],
                       z = c(-1,-1,1,1,-1,-1,1,1) / 2 * aspect[2])

        ## these are box boundaries
        ##pre <- c(1,2,3,4,2,3,4,1,5,6,7,8)
        ##nxt <- c(2,3,4,1,6,7,8,5,6,7,8,5)  old -- as backup

        pre <- c(1,2,4,1,2,3,4,1,5,6,8,5)
        nxt <- c(2,3,3,4,6,7,8,5,6,7,7,8)

        ## SCALES : very beta

        labs <- rbind(x = c(0, corners$x[pre[scpos$y]], corners$x[pre[scpos$z]]),
                      y = c(corners$y[pre[scpos$x]], 0, corners$y[pre[scpos$z]]),
                      z = c(corners$z[pre[scpos$x]], corners$z[pre[scpos$y]], 0))

        labs[,1] <- labs[,1] * (1 + scales.3d$x.scales$distance/3)
        labs[,2] <- labs[,2] * (1 + scales.3d$y.scales$distance/3)
        labs[,3] <- labs[,3] * (1 + scales.3d$z.scales$distance/3)

        axes <- rbind(x = 
                      c(proportion * corners$x[c(pre[scpos$x], nxt[scpos$x])],
                        corners$x[c(pre[scpos$y], nxt[scpos$y])],
                        corners$x[c(pre[scpos$z], nxt[scpos$z])]),
                      y = 
                      c(corners$y[c(pre[scpos$x], nxt[scpos$x])],
                        proportion * corners$y[c(pre[scpos$y], nxt[scpos$y])],
                        corners$y[c(pre[scpos$z], nxt[scpos$z])]),
                      z = 
                      c(corners$z[c(pre[scpos$x], nxt[scpos$x])],
                        corners$z[c(pre[scpos$y], nxt[scpos$y])],
                        proportion * corners$z[c(pre[scpos$z], nxt[scpos$z])]))
            
        axes[,1:2] <- axes[,1:2] * (1 + scales.3d$x.scales$distance/10)
        axes[,3:4] <- axes[,3:4] * (1 + scales.3d$y.scales$distance/10)
        axes[,5:6] <- axes[,5:6] * (1 + scales.3d$z.scales$distance/10)

        x.at <-
            if (is.logical(scales.3d$x.scales$at))
                lpretty(xlim, scales.3d$x.scales$tick.number)
            else scales.3d$x.scales$at
        y.at <- 
            if (is.logical(scales.3d$y.scales$at))
                lpretty(ylim, scales.3d$y.scales$tick.number)
            else scales.3d$y.scales$at
        z.at <- 
            if (is.logical(scales.3d$z.scales$at))
                lpretty(zlim, scales.3d$z.scales$tick.number)
            else scales.3d$z.scales$at
        x.at <- x.at[x.at >= xlim[1] & x.at < xlim[2]]
        y.at <- y.at[y.at >= ylim[1] & y.at < ylim[2]]
        z.at <- z.at[z.at >= zlim[1] & z.at < zlim[2]]
        x.at.lab <-
            if (is.logical(scales.3d$x.scales$labels))
                as.character(x.at)
            else as.character(scales.3d$x.scales$labels)
        y.at.lab <-
            if (is.logical(scales.3d$y.scales$labels))
                as.character(y.at)
            else as.character(scales.3d$y.scales$labels)
        z.at.lab <-
            if (is.logical(scales.3d$z.scales$labels))
                as.character(z.at)
            else as.character(scales.3d$z.scales$labels)

    

        ## box ranges and lengths
        cmin <- lapply(corners, min)
        cmax <- lapply(corners, max)
        clen <- lapply(corners, function(x) diff(range(x)))

        ##if (!missing(R.mat))
        ##  {;} #do something

        ## transformed data
        tdata <- rbind(x = cmin$x + clen$x * (x-xlim[1])/diff(xlim),
                       y = cmin$y + clen$y * (y-ylim[1])/diff(ylim),
                       z = cmin$z + clen$z * (z-zlim[1])/diff(zlim))

        taxes <- rot.mat %*% axes
        x.at <- cmin$x + clen$x * (x.at-xlim[1])/diff(xlim)
        y.at <- cmin$y + clen$y * (y.at-ylim[1])/diff(ylim)
        z.at <- cmin$z + clen$z * (z.at-zlim[1])/diff(zlim)
        at.len <- length(x.at)
        x.at <- rbind(x = x.at,
                      y = rep(corners$y[pre[scpos$x]], at.len),
                      z = rep(corners$z[pre[scpos$x]], at.len))
        at.len <- length(y.at)
        y.at <- rbind(x = rep(corners$x[pre[scpos$y]], at.len),
                      y = y.at,
                      z = rep(corners$z[pre[scpos$y]], at.len))
        at.len <- length(z.at)
        z.at <- rbind(x = rep(corners$x[pre[scpos$z]], at.len),
                      y = rep(corners$y[pre[scpos$z]], at.len),
                      z = z.at)

        x.at.end <- x.at + scales.3d$x.scales$tck * .05 * labs[,1]
        y.at.end <- y.at + scales.3d$y.scales$tck * .05 * labs[,2]
        z.at.end <- z.at + scales.3d$z.scales$tck * .05 * labs[,3]

        x.labs <- x.at + 2 * scales.3d$x.scales$tck * .05 * labs[,1]
        y.labs <- y.at + 2 * scales.3d$y.scales$tck * .05 * labs[,2]
        z.labs <- z.at + 2 * scales.3d$z.scales$tck * .05 * labs[,3]

        x.at <- rot.mat %*% x.at
        x.labs <- rot.mat %*% x.labs
        x.at.end <- rot.mat %*% x.at.end
        y.at <- rot.mat %*% y.at
        y.labs <- rot.mat %*% y.labs
        y.at.end <- rot.mat %*% y.at.end
        z.at <- rot.mat %*% z.at
        z.labs <- rot.mat %*% z.labs
        z.at.end <- rot.mat %*% z.at.end
        
        tdata <- rot.mat %*% tdata
        corners <- rot.mat %*% t(as.matrix(corners))
        tlabs <- rot.mat %*% labs


        zback <- min(corners[3,])
        zfront <- max(corners[3,])
        za <- (zfront * (1-distance) - zback) / (zfront - zback)
        zb <- distance / (zfront - zback)

        tdata[1,] <- (za + zb * tdata[3,]) * tdata[1,]
        tdata[2,] <- (za + zb * tdata[3,]) * tdata[2,]

        corners[1,] <- (za + zb * corners[3,]) * corners[1,]
        corners[2,] <- (za + zb * corners[3,]) * corners[2,]

        taxes[1,] <- (za + zb * taxes[3,]) * taxes[1,]
        taxes[2,] <- (za + zb * taxes[3,]) * taxes[2,]
        
        x.at[1,] <- (za + zb * x.at[3,]) * x.at[1,]
        x.at[2,] <- (za + zb * x.at[3,]) * x.at[2,]
        x.labs[1,] <- (za + zb * x.labs[3,]) * x.labs[1,]
        x.labs[2,] <- (za + zb * x.labs[3,]) * x.labs[2,]
        x.at.end[1,] <- (za + zb * x.at.end[3,]) * x.at.end[1,]
        x.at.end[2,] <- (za + zb * x.at.end[3,]) * x.at.end[2,]

        y.at[1,] <- (za + zb * y.at[3,]) * y.at[1,]
        y.at[2,] <- (za + zb * y.at[3,]) * y.at[2,]
        y.labs[1,] <- (za + zb * y.labs[3,]) * y.labs[1,]
        y.labs[2,] <- (za + zb * y.labs[3,]) * y.labs[2,]
        y.at.end[1,] <- (za + zb * y.at.end[3,]) * y.at.end[1,]
        y.at.end[2,] <- (za + zb * y.at.end[3,]) * y.at.end[2,]

        z.at[1,] <- (za + zb * z.at[3,]) * z.at[1,]
        z.at[2,] <- (za + zb * z.at[3,]) * z.at[2,]
        z.labs[1,] <- (za + zb * z.labs[3,]) * z.labs[1,]
        z.labs[2,] <- (za + zb * z.labs[3,]) * z.labs[2,]
        z.at.end[1,] <- (za + zb * z.at.end[3,]) * z.at.end[1,]
        z.at.end[2,] <- (za + zb * z.at.end[3,]) * z.at.end[2,]

        
        tlabs[1,] <- (za + zb * tlabs[3,]) * tlabs[1,]
        tlabs[2,] <- (za + zb * tlabs[3,]) * tlabs[2,]


        farthest <- 1
        farval <- corners[3,1]

        for (i in 2:8)
            if (corners[3,i] < farval) {
                farthest <- i
                farval <- corners[3,i]
            }

        mark <- rep(TRUE, 12)
        for (j in 1:12)
            if (pre[j]==farthest || nxt[j]==farthest)
                mark[j] <- FALSE

        lsegments(corners[1, pre[!mark]],
                  corners[2, pre[!mark]],
                  corners[1, nxt[!mark]],
                  corners[2, nxt[!mark]],
                  col = par.box.final$col,
                  lwd = par.box.final$lwd,
                  lty = 2)

        
        if (wireframe) {
            ## This is where the wireframe is actually drawn
            if (TRUE) {
                xx <- tdata[1,]
                yy <- tdata[2,]
                zz <- tdata[3,]
                ord <- order(zz[id0])
                zcol <- zcol[id0][ord]
                px <- cbind(xx[id0][ord], xx[id1][ord], xx[id2][ord], xx[id3][ord])
                py <- cbind(yy[id0][ord], yy[id1][ord], yy[id2][ord], yy[id3][ord])
                for (i in seq(along = ord))
                    grid.polygon(x = px[i,], y = py[i,], default.units = "native",
                                 gp = gpar(fill = col.regions[zcol[i]], col = "black"))
            }
#             else {
#                 xx <- tdata[1,]
#                 yy <- tdata[2,]
#                 zz <- tdata[3,]
#                 ord <- order(zz[id0])

#                 dumx <- unit(c(0,.5,.5,0 ), "native")
#                 dumy <- unit(c(0,0,.5,.5), "native")

#                 grid.Call.graphics("D_quadrilateral",
#                                    xx[id0][ord],
#                                    xx[id1][ord],
#                                    xx[id2][ord],
#                                    xx[id3][ord],
#                                    yy[id0][ord],
#                                    yy[id1][ord],
#                                    yy[id2][ord],
#                                    yy[id3][ord],
#                                    dumx, dumy,
#                                    col.regions[zcol[id0][ord]],
#                                    PACKAGE="grid")
#             }
        }
        else subpanel(x=tdata[1,], y=tdata[2,], subscripts = subscripts, ...)
            
        lsegments(corners[1, pre[mark]],
                  corners[2, pre[mark]],
                  corners[1, nxt[mark]],
                  corners[2, nxt[mark]],
                  col = par.box.final$col,
                  lty = par.box.final$lty,
                  lwd = par.box.final$lwd)

        ## Next part for axes : beta
        
        if (scales.3d$x.scales$draw) {
            if (scales.3d$x.scales$arrows) {
                larrows(x0 = taxes[1, 1], y0 = taxes[2, 1],
                        x1 = taxes[1, 2], y1 = taxes[2, 2],
                        lty = scales.3d$x.scales$lty,
                        lwd = scales.3d$x.scales$lwd,
                        col = scales.3d$x.scales$col)
            }
            else {
                lsegments(x0 = x.at[1,], y0 = x.at[2,], x1 = x.at.end[1,], y1 = x.at.end[2,],
                          lty = scales.3d$x.scales$lty,
                          col = scales.3d$x.scales$col,
                          lwd = scales.3d$x.scales$lwd)
                ltext(x.at.lab, x = x.labs[1,], y = x.labs[2,],
                      cex = scales.3d$x.scales$cex,
                      font = scales.3d$x.scales$font,
                      col = scales.3d$x.scales$col)
            }
        }

        if (scales.3d$y.scales$draw) {
            if (scales.3d$y.scales$arrows) {
                larrows(x0 = taxes[1, 3], y0 = taxes[2, 3],
                        x1 = taxes[1, 4], y1 = taxes[2, 4],
                        lty = scales.3d$y.scales$lty,
                        lwd = scales.3d$y.scales$lwd,
                        col = scales.3d$y.scales$col)
            }
            else {
                lsegments(x0 = y.at[1,], y0 = y.at[2,], x1 = y.at.end[1,], y1 = y.at.end[2,],
                          lty = scales.3d$y.scales$lty,
                          col = scales.3d$y.scales$col,
                          lwd = scales.3d$y.scales$lwd)
                ltext(y.at.lab, x = y.labs[1,], y = y.labs[2,],
                      cex = scales.3d$y.scales$cex,
                      font = scales.3d$y.scales$font,
                      col = scales.3d$y.scales$col)
            }
        }
        if (scales.3d$z.scales$draw) {
            if (scales.3d$z.scales$arrows) {
                larrows(x0 = taxes[1, 5], y0 = taxes[2, 5],
                        x1 = taxes[1, 6], y1 = taxes[2, 6],
                        lty = scales.3d$z.scales$lty,
                        lwd = scales.3d$z.scales$lwd,
                        col = scales.3d$z.scales$col)
            }
            else {
                lsegments(x0 = z.at[1,], y0 = z.at[2,], x1 = z.at.end[1,], y1 = z.at.end[2,],
                          lty = scales.3d$z.scales$lty,
                          col = scales.3d$x.scales$col,
                          lwd = scales.3d$z.scales$lwd)
                ltext(z.at.lab, x = z.labs[1,], y = z.labs[2,],
                      cex = scales.3d$z.scales$cex,
                      font = scales.3d$z.scales$font,
                      col = scales.3d$z.scales$col)
            }
        }


        if (!is.null(xlab)) ltext(xlab$lab, x = tlabs[1, 1], y = tlabs[2, 1],
                                  cex = xlab$cex, rot = xlab$rot,
                                  font = xlab$font, col = xlab$col)
        
        if (!is.null(ylab)) ltext(ylab$lab, x = tlabs[1, 2], y = tlabs[2, 2],
                                  cex = ylab$cex, rot = ylab$rot,
                                  font = ylab$font, col = ylab$col)
                                  
        if (!is.null(zlab)) ltext(zlab$lab, x = tlabs[1, 3], y = tlabs[2, 3],
                                  cex = zlab$cex, rot = zlab$rot,
                                  font = zlab$font, col = zlab$col)
    }
}








panel.wireframe <- function(...)
    panel.cloud(..., wireframe = TRUE)





wireframe <-
    function(formula,
             data = parent.frame(),
             panel = "panel.wireframe",
             prepanel = NULL,
             strip = TRUE,
             groups = NULL,
             cuts = 70,
             pretty = FALSE,
             col.regions = trellis.par.get("regions")$col,
             drape = FALSE,
             colorkey = any(drape),
             ...,
             subset = TRUE)
{
    warning("wireframe can be EXTREMELY slow")
    ## m <- match.call(expand.dots = FALSE)
    dots <- list(...)
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    do.call("cloud",
            c(list(formula = formula, data = data,
                   groups = groups, subset = subset,
                   panel = panel, prepanel = prepanel, strip = strip,
                   cuts = cuts, 
                   pretty = pretty,
                   col.regions = col.regions,
                   drape = drape,
                   colorkey = colorkey),
              dots))
}























cloud <-
    function(formula,
             data = parent.frame(),
             aspect = c(1,1),
             layout = NULL,
             panel = "panel.cloud",
             subpanel = "panel.xyplot",
             prepanel = NULL,
             scales = NULL,
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim = range(x),
             ylab,
             ylim = range(y),
             zlab,
             zlim = range(z),
             distance = .2,
             par.box,
             perspective = TRUE,
             R.mat = diag(4),
             screen = list(z = 40, x = -60),
             zoom = .9,
             at,
             pretty = FALSE,
             drape = FALSE,
             colorkey = any(drape),
             col.regions, cuts = 1,
             ...,
             subscripts = TRUE,
             subset = TRUE)
{

    ##dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, z etc. and do some preprocessing


    formula <- eval(substitute(formula), data, parent.frame())
    form <-
        if (inherits(formula, "formula"))
            latticeParseFormula(formula, data, dim = 3)
        else {
            if (!is.matrix(formula)) stop("invalid formula")
            else {
                tmp <- expand.grid(1:nrow(formula), 1:ncol(formula))
                list(left = as.vector(formula),
                     right.x = tmp[[1]],
                     right.y = tmp[[2]],
                     condition = NULL,
                     left.name = "",
                     right.x.name = "", right.y.name = "")
            }
        }


    cond <- form$condition
    number.of.cond <- length(cond)
    z <- form$left
    x <- form$right.x
    y <- form$right.y

    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    y <- y[subset, drop = TRUE]
    z <- z[subset, drop = TRUE]
    subscr <- subscr[subset, drop = TRUE]
    
    if (missing(xlab)) xlab <- form$right.x.name
    if (missing(ylab)) ylab <- form$right.y.name
    if (missing(zlab)) zlab <- form$left.name

    if(!(is.numeric(x) && is.numeric(y) && is.numeric(z)))
        warning("x, y and z should be numeric")
    x <- as.numeric(x)
    y <- as.numeric(y)
    z <- as.numeric(z)

    zrng <- extend.limits(range(z[!is.na(z)]))
    if (missing(at))
        at <-
            if (pretty) lpretty(zrng, cuts)
            else seq(zrng[1], zrng[2], length = cuts+2)
    

    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = 1,
                          strip = strip,
                          panel = panel,
                          xlab = NULL,
                          ylab = NULL), dots))
                          
    ##-----------------------------------------------------------
    ## xlab, ylab, zlab have special meaning in cloud / wireframe

    if (!is.null(xlab)) {
        text <- trellis.par.get("par.xlab.text")
        if (is.null(text)) text <- list(cex = 1, col = "black", font = 1)
        xlab <- list(label = xlab[[1]], col = text$col, rot = 0,
                     cex = text$cex, font = text$font)
        if (is.list(xlab)) xlab[names(xlab)] <- xlab
    }
    if (!is.null(ylab)) {
        text <- trellis.par.get("par.ylab.text")
        if (is.null(text)) text <- list(cex = 1, col = "black", font = 1)
        ylab <- list(label = ylab[[1]], col = text$col,  rot = 0,
                         cex = text$cex, font = text$font)
        if (is.list(ylab)) ylab[names(ylab)] <- ylab
    }
    if (!is.null(zlab)) {
        text <- trellis.par.get("par.zlab.text")
        if (is.null(text)) text <- list(cex = 1, col = "black", font = 1)
        zlab <- list(label = zlab[[1]], col = text$col, rot = 0,
                         cex = text$cex, font = text$font)
        if (is.list(zlab)) zlab[names(zlab)] <- zlab
    }
    ##-----------------------------------------------------------------





    
    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(xlab) && !is.character(xlab$label))
        xlab$label <- form$right.x.name
    if (is.list(ylab) && !is.character(ylab$label))
        ylab$label <- form$right.y.name
    if (is.list(zlab) && !is.character(zlab$label))
        zlab$label <- form$left.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    foo <- c(foo,
             do.call("construct.scales", list(draw=FALSE)))

    ## scales has to be interpreted differently. Nothing needs to be
    ## done for the ususal scales, but need a scales for panel.cloud
    ## Splus probably doesn't allow x-y-z-specific scales, but I see
    ## no reason not to allow that (will not allow limits, though)


    scales.default <-
        list(cex = .8, col = "black", lty = 1,
             lwd = 1, tck = 1, distance = c(1, 1, 1),
             arrows = TRUE)
    if (!is.null(scales)) scales.default[names(scales)] <- scales
    scales.3d <- do.call("construct.3d.scales", scales.default)

    ## Step 3: Decide if limits were specified in call:
    ## Here, always FALSE (in the 2d panel sense)
    have.xlim <- FALSE
    have.ylim <- FALSE

    ## Step 4: Decide if log scales are being used: !!!

    have.xlog <- !is.logical(scales.3d$x.scales$log) || scales.3d$x.scales$log
    have.ylog <- !is.logical(scales.3d$y.scales$log) || scales.3d$y.scales$log
    have.zlog <- !is.logical(scales.3d$z.scales$log) || scales.3d$z.scales$log
    if (have.xlog) {
        xlog <- scales.3d$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        ylog <- scales.3d$y.scales$log
        ybase <-
            if (is.logical(ylog)) 10
            else if (is.numeric(ylog)) ylog
            else if (ylog == "e") exp(1)

        y <- log(y, ybase)
        ylim <- log(ylim, ybase)
    }
    if (have.zlog) {
        zlog <- scales.3d$z.scales$log
        zbase <-
            if (is.logical(zlog)) 10
            else if (is.numeric(zlog)) zlog
            else if (zlog == "e") exp(1)

        z <- log(z, zbase)
        zlim <- log(zlim, zbase)
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)|is.na(y)|is.na(z)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    ## FIXME: R.mat ignored
    ## calculate rotation matrix:
    rot.mat <- diag(3)
    screen.names <- names(screen)
    screen <- lapply(screen, "*", pi/180)

    for(i in seq(along=screen.names)) {
        th <- screen[[i]]
        cth <- cos(th)
        sth <- sin(th)
        tmp.mat <- 
            (if (screen.names[i]=="x")
             matrix(c(1, 0, 0, 0, cth, sth, 0, -sth, cth), 3, 3)
            else if (screen.names[i]=="y")
             matrix(c(cth, 0, -sth, 0, 1, 0, sth, 0, cth), 3, 3)
            else if (screen.names[i]=="z")
             matrix(c(cth, sth, 0, -sth, cth, 0, 0, 0, 1), 3, 3))
        rot.mat <- tmp.mat %*% rot.mat
    }


    if (drape) {
        ## region
        numcol <- length(at)-1
        numcol.r <- length(col.regions)

        col.regions <-
            if (numcol.r <= numcol)
                rep(col.regions, length = numcol)
            else col.regions[floor(1+(1:numcol-1)*(numcol.r-1)/(numcol-1))]
    
        if (is.logical(colorkey)) {
            if (colorkey) foo$colorkey <-
                list(space = "right", col = col.regions,
                     at = at, tick.number = 7)
        }
        else if (is.list(colorkey)) {
            foo$colorkey <- colorkey
            if (is.null(foo$colorkey$col)) foo$colorkey$col <- col.regions
            if (is.null(foo$colorkey$at)) foo$colorkey$at <- at
        }
        zcol <- numeric(length(z))
        for (i in seq(along=col.regions))
            zcol[!id.na & z>=at[i] & z<at[i+1]] <- i
    }
    else {
        col.regions <- trellis.par.get("background")$col
        zcol <- 1
    }









    ## maybe *lim = NULL with relation = "free" ? 
    foo$panel.args.common <-
        c(list(x=x, y=y, z=z, rot.mat = rot.mat, zoom = zoom,
               subpanel = subpanel,
               xlim = xlim, ylim = ylim, zlim = zlim,
               xlab = xlab, ylab = ylab, zlab = zlab,
               aspect = aspect,
               distance = distance,
               scales.3d = scales.3d,
               zcol = zcol, col.regions = col.regions),
          dots)

    if (!is.null(groups)) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number

    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    foo$panel.args[[panel.number]] <- 
                        list(subscripts = subscr[id])

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.cloud,
                               prepanel = prepanel,
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = 1,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}










### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA





cupdate <- function(index, maxim)
{
    ##  This function is used to handle arbitrary number of
    ## conditioning variables : every time it is called, it
    ## increments the "current" level of the conditioning
    ## variables suitably, i.e., it tries to increment the
    ## level of the 1st conditining variable (the one which
    ## varies fastest along panel order) and if it happens
    ## to be at its maximum (last) value, it sets it to the
    ## first value AND increments the "current" level of the
    ## 2nd (next) conditioning variable recursively.

    ## This is an internal function, not to be documented
    ## for the high level user.
    
    if(length(index)!=length(maxim)||length(maxim)<=0)
        stop("Inappropriate arguments")
    index[1] <- index[1] + 1
    if(index[1]>maxim[1] && length(maxim)>1)
        c(1,cupdate(index[-1],maxim[-1]))
    else index
}




latticeParseFormula <-
    function(model, data, dimension = 2)
{
    parseCond <-
        function(model)
        {
            model <- eval(parse(text = paste("~", deparse(model))))[[2]]
            model.vars <- list()
            while (length(model) == 3 && (model[[1]] == as.name("*")
                         || model[[1]] == as.name("+"))) {
                model.vars <- c(model.vars, model[[3]])
                model <- model[[2]]
            }
            rev(c(model.vars, model))
        }
    if (!inherits(model, "formula"))
        stop("model must be a formula object")

    ans <- if (dimension == 2) {
        list(left = NULL, right = NULL, condition = NULL,
             left.name = character(0), right.name = character(0))
    }
    else if (dimension == 3) {
        list(left = NULL, right.x = NULL, right.y = NULL, condition = NULL,
             left.name = character(0), right.x.name = character(0),
             right.y.name = character(0))
    }
    else stop(paste("invalid dimension : ", dimension))

    if (length(model) == 3) {
        ans$left <- eval(model[[2]], data)
        ans$left.name <- deparse(model[[2]])
    }
    model <- model[[length(model)]]
    if (length(model) == 3 && model[[1]] == as.name("|")) {
        model.vars <- parseCond(model[[3]])
        ans$condition <- vector("list", length(model.vars))
        names(ans$condition) <- sapply(model.vars, deparse)
        for (i in seq(along = model.vars))
            ans$condition[[i]] <- eval(model.vars[[i]], data)
        model <- model[[2]]
    }
    if (dimension == 2) {
        ans$right <- eval(model, data)
        ans$right.name <- deparse(model)
    }
    else if (dimension == 3 && length(model) == 3 &&
             (model[[1]] == "*" || model[[1]] == "+")) {
        ans$right.x <- eval(model[[2]], data)
        ans$right.y <- eval(model[[3]], data)
        ans$right.x.name <- deparse(model[[2]])
        ans$right.y.name <- deparse(model[[3]])
    }
    else stop("invalid model")

    ans
}





banking <- function(dx, dy)
{
    if (is.na(dx)) NA
    else {
        if (is.list(dx)) {
            dy <- dx[[2]]
            dx <- dx[[1]]
        }
        if (length(dx)!=length(dy)) stop("Non matching lengths")
        id <- dx!=0 & dy!=0
        if (any(id)) {
            r  <- abs(dx[id]/dy[id])
            median(r)
        }
        else 1
    }
}





extend.limits <-
    function(lim, length=1, prop = 0.07) {
        if (!is.numeric(lim)) NA
        else if(length(lim)==2) {
            if (lim[1]>lim[2]) stop("Improper value of limit")
            if(!missing(length) && !missing(prop))
                stop("length and prop cannot both be specified")
            if(length <= 0) stop("length must be positive")
            if(!missing(length)) prop <- (length-diff(lim))/(2*diff(lim))
            if(lim[1]==lim[2]) lim + 0.5*c(-length,length)
            else {
                d <- diff(lim)
                lim + prop*d*c(-1,1)
            }
        }
        else {
            print(lim)
            stop("improper length of lim in extend.limits")
        }
    }






construct.scales <-
    function(draw = TRUE,
             tck = 1,
             tick.number = 5,
             cex = 1,
             rot = FALSE,
             at = FALSE,
             labels = FALSE,
             col = FALSE,
             log = FALSE,
             font = FALSE,
             alternating = TRUE,
             relation = "same",
             x = NULL,
             y = NULL,
             ...)
{
    xfoo <- list(draw = draw, tck = tck,
                 tick.number = tick.number,
                 cex = cex, rot = rot, font = font,
                 at = at, labels = labels,
                 col = col, log = log,
                 alternating = alternating,
                 relation = relation)
    yfoo <- xfoo
    if (!is.null(x)) {
        if (is.character(x)) x <- list(relation = x)
        xfoo[names(x)] <- x
    }
    if (is.logical(xfoo$alternating))
        xfoo$alternating <-
            if (xfoo$alternating) c(1,2)
            else 1
    if (!is.null(y)) {
        if (is.character(y)) y <- list(relation = y)
        yfoo[names(y)] <- y
    }
    if (is.logical(yfoo$alternating))
        yfoo$alternating <-
            if (yfoo$alternating) c(1,2)
            else 1
    list(x.scales = xfoo, y.scales = yfoo)
}





construct.3d.scales <-
    function(draw = TRUE,
             tck = 1,
             lty = 1, lwd = 1,
             distance = c(1,1,1),
             tick.number = 5,
             cex = 1,
             rot = FALSE,
             at = FALSE,
             labels = FALSE,
             col = FALSE,
             log = FALSE,
             font = FALSE,
             arrows = TRUE,
             relation = "same",
             x = NULL,
             y = NULL,
             z = NULL,
             ...)
{
    xfoo <- list(draw = draw, tck = tck,
                 lty = 1, lwd = 1,
                 tick.number = tick.number,
                 cex = cex, rot = rot, font = font,
                 at = at, labels = labels,
                 col = col, log = log, arrows = arrows,
                 relation = relation)
    yfoo <- xfoo
    zfoo <- xfoo
    xfoo$distance <- distance[1]
    yfoo$distance <- distance[2]
    zfoo$distance <- distance[3]
    if (!is.null(x)) {
        if (is.character(x)) x <- list(relation = x)
        xfoo[names(x)] <- x
    }
    if (!is.null(y)) {
        if (is.character(y)) y <- list(relation = y)
        yfoo[names(y)] <- y
    }
    if (!is.null(z)) {
        if (is.character(z)) z <- list(relation = z)
        zfoo[names(z)] <- z
    }
    list(x.scales = xfoo, y.scales = yfoo, z.scales = zfoo)
}




limits.and.aspect <-
    function(prepanel.default.function,
             prepanel = NULL,
             have.xlim = FALSE, xlim = NULL,
             have.ylim = FALSE, ylim = NULL,
             x.relation, y.relation,
             panel.args.common = list(),
             panel.args = list(),
             aspect,
             nplots,
             ...)  ## extra arguments for prepanel (for qqmathline)
{

    if (nplots<1) stop("need at least one panel")
    x.limits <- as.list(1:nplots)
    y.limits <- as.list(1:nplots)
    dxdy <- as.list(1:nplots)


    for (count in 1:nplots)
    {
        if (is.list(panel.args[[count]])) {

            pargs <- c(panel.args.common, panel.args[[count]], list(...))
            tem <- do.call("prepanel.default.function", pargs)
            if (is.function(prepanel)) {
                prenames <- names(formals(prepanel))
                if (!("..." %in% prenames)) pargs <- pargs[prenames]
                pretem <- do.call("prepanel", pargs)
                tem[names(pretem)] <- pretem
            }
            x.limits[[count]] <- tem$xlim
            y.limits[[count]] <- tem$ylim
            dxdy[[count]] <- list(tem$dx, tem$dy)

        }
        else {
            x.limits[[count]] <- c(NA, NA)
            y.limits[[count]] <- c(NA, NA)
            dxdy[[count]] <- NA
        } 

    }

    if (x.relation == "same")
        x.limits <-
            if (have.xlim) xlim
            else extend.limits(range(unlist(x.limits), na.rm = TRUE))
    else if (x.relation == "sliced") {
        x.slicelen <- 1.14 * max(unlist(lapply(x.limits, diff)), na.rm = TRUE)
        x.limits <- lapply(x.limits, extend.limits, length = x.slicelen)
    }
    else if (x.relation == "free")
        x.limits <- lapply(x.limits, extend.limits)

    if (y.relation == "same") {
        y.limits <-
            if (have.ylim) ylim
            else extend.limits(range(unlist(y.limits), na.rm = TRUE))
    }
    else if (y.relation == "sliced") {
        y.slicelen <- 1.14 * max(unlist(lapply(y.limits, diff)), na.rm = TRUE)
        y.limits <- lapply(y.limits, extend.limits, length = y.slicelen)
    }
    else if (y.relation == "free")
        y.limits <- lapply(y.limits, extend.limits)

    if (is.character(aspect))
        if (aspect == "xy") {
            aspect <- median(unlist(lapply(dxdy, banking)), na.rm = TRUE)
            if (y.relation == "free")
                warning("aspect=xy when y-relation=free is not sensible")
            else aspect <- aspect *
                if (y.relation == "sliced") y.slicelen else diff(y.limits)
            if (x.relation == "free")
                warning("aspect=xy when x-relation=free is not sensible")
            else aspect <- aspect /
                if (x.relation == "sliced") x.slicelen else diff(x.limits)
        }
    else aspect <- 1

    list(x.limits = x.limits,
         y.limits = y.limits,
         aspect.ratio = aspect)
}






trellis.skeleton <-
    function(as.table = FALSE,
             aspect = "fill",
             between = list(x=0, y=0),
             key = NULL,
             page = NULL,
             main = NULL,
             sub = NULL,
             par.strip.text = list(),
             skip = FALSE,
             strip = strip.default,
             xlab = NULL,
             ylab = NULL,
             panel = panel,
             ...)
{
    foo <- list(as.table = as.table,
                aspect.fill = aspect=="fill",
                key = key,
                panel = panel, 
                page = page,
                skip = skip,
                strip = if (is.logical(strip) && strip) "strip.default"
                else strip,
                x.between = 0,
                y.between = 0,
                par.strip.text = trellis.par.get("add.text"))
    if (is.null(foo$par.strip.text)) foo$par.strip.text = list(col = "black", cex = 1, font = 1)
    foo$par.strip.text$lines <- 1
    
    if (!is.null(between$x)) foo$x.between <- between$x
    if (!is.null(between$y)) foo$y.between <- between$y

    foo$par.strip.text[names(par.strip.text)] <- par.strip.text

    if (!is.null(main)) {
        text <- trellis.par.get("par.main.text")
        if (is.null(text)) text <- list(cex = 1.2, col = "black", font = 2)
        foo$main <- list(label = main[[1]], col = text$col, cex = text$cex, font = text$font)
        if (is.list(main)) foo$main[names(main)] <- main
    }
    if (!is.null(sub)) {
        text <- trellis.par.get("par.sub.text")
        if (is.null(text)) text <- list(cex = 1, col = "black", font = 2)
        foo$sub <- list(label = sub[[1]], col = text$col, cex = text$cex, font = text$font)
        if (is.list(sub)) foo$sub[names(sub)] <- sub
    }
    if (!is.null(xlab)) {
        text <- trellis.par.get("par.xlab.text")
        if (is.null(text)) text <- list(cex = 1, col = "black", font = 1)
        foo$xlab <- list(label = xlab[[1]], col = text$col, cex = text$cex, font = text$font)
        if (is.list(xlab)) foo$xlab[names(xlab)] <- xlab
    }
    if (!is.null(ylab)) {
        text <- trellis.par.get("par.ylab.text")
        if (is.null(text)) text <- list(cex = 1.2, col = "black", font = 2)
        foo$ylab <- list(label = ylab[[1]], col = text$col, cex = text$cex, font = text$font)
        if (is.list(ylab)) foo$ylab[names(ylab)] <- ylab
    }
list(foo = foo, dots = list(...))
}







compute.layout <-
    function(layout, cond.max.level, skip = FALSE)
{
    number.of.cond <- length(cond.max.level)
    nplots <- prod(cond.max.level)
    
    if (!is.numeric(layout)) {
        layout <- c(0,1,1)
        if (number.of.cond==1) layout[2] <- nplots
        else {
            layout[1] <- cond.max.level[1]
            layout[2] <- cond.max.level[2]
        }
        skip <- rep(skip, length = max(layout[1] * layout[2], layout[2]))
        plots.per.page <- length(skip) - length(skip[skip])
        layout[3] <- ceiling(nplots/plots.per.page) # + 1
    }
    else if (length(layout)==1)
        stop("layout must have at least 2 elements")
    else if (length(layout)==2)
    {
        if(all(layout<1))
            stop("at least one element of layout must be positive")
        else if (layout[2]==0) stop("inadmissible value of layout")
        
        skip <- rep(skip, length = max(layout[1] * layout[2], layout[2]))
        plots.per.page <- length(skip) - length(skip[skip])
        layout[3] <- ceiling(nplots/plots.per.page) # + 1 
    }
    else if (length(layout)==3) {
        if(layout[1]<0||layout[2]<1||layout[3]<1)
            stop("invalid value for layout")
    }
    layout
}





### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA









prepanel.default.densityplot <-
    function(x,
             darg,
             ...)
{
    if (length(x)<1)
        list(xlim = NA,
             ylim = NA,
             dx = NA,
             dy = NA)
    else
    {
        h <- do.call("density", c(list(x=x), darg))
        list(xlim = range(h$x),
             ylim = range(h$y),
             dx = diff(h$x), dy = diff(h$y))
    }
}




panel.densityplot <-
    function(x,
             darg = list(n = 30),
             plot.points = TRUE,
             ref = FALSE,
             cex = 0.5,
             col = plot.line$col,
             col.line,
             ...)
{
    if (ref) {
        reference.line <- trellis.par.get("reference.line")
        panel.abline(h=0,
                     col = reference.line$col,
                     lty = reference.line$lty,
                     lwd = reference.line$lwd)
    }
    if (length(x)>1) {
        plot.line <- trellis.par.get("plot.line")
        if (missing(col.line)) col.line <- col
        h <- do.call("density", c(list(x=x), darg))
        lim <- current.viewport()$xscale
        id <- (h$x>=lim[1] & h$x<=lim[2])
        llines(x = h$x[id], y = h$y[id], col = col.line, ...)
        if (plot.points) panel.xyplot(x = x, y = rep(0, length(x)), cex = cex, col = col, ...) 
    }
}





densityplot <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.densityplot",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             bw = NULL,
             adjust = NULL,
             kernel = NULL,
             window = NULL,
             width = NULL,
             give.Rkern = FALSE,
             n = 50,
             from = NULL,
             to = NULL,
             cut = NULL,
             na.rm = NULL,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ## dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## darg is a list that gives arguments to density()
    darg <- list()
    darg$bw <- bw
    darg$adjust <- adjust
    darg$kernel <- kernel
    darg$window <- window
    darg$width <- width
    darg$give.Rkern <- give.Rkern
    darg$n <- n
    darg$from <- from
    darg$to <- to
    darg$cut <- cut
    darg$na.rm <- na.rm
    
    ## Step 1: Evaluate x, y, etc. and do some preprocessing
    
    formname <- deparse(substitute(formula))
    formula <- eval(substitute(formula), data, parent.frame())

    form <-
        if (inherits(formula, "formula"))
            latticeParseFormula(formula, data)
        else {
            if (!is.numeric(formula)) stop("invalid formula")
            else {
                list(left = NULL,
                     right = formula,
                     condition = NULL,
                     left.name = "",
                     right.name = formname)
            }
        }

    ##form <- latticeParseFormula(formula, data)



    cond <- form$condition
    number.of.cond <- length(cond)
    x <- form$right
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    if(subscripts) subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    if (subscripts) subscr <- subscr[subset, drop = TRUE]
    
    if (missing(xlab)) xlab <- form$right.name
    if (missing(ylab)) ylab <- "Density"

    if (!is.numeric(x))
        warning("x should be numeric")
    x <- as.numeric(x)

    ## create a skeleton trellis object with the
    ## less complicated components:
    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- form$right.name
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <- "Density"

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo,
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:
    
    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)
        
        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        warning("Can't have log Y-scale")
        have.ylog <- FALSE
        foo$y.scales$log <- FALSE
    }

    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    foo$panel.args.common <- c(dots, list(darg = darg))
    if (subscripts) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    foo$panel.args[[panel.number]] <-
                        list(x = x[id])
                    if (subscripts)
                        foo$panel.args[[panel.number]]$subscripts <-
                            subscr[id]

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.densityplot,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}


### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA









prepanel.default.histogram <-
    function(x,
             breaks = NULL,
             equal.widths = TRUE,
             type = "density",
             ...)
{
    if (length(x)<1)
        list(xlim = NA,
             ylim = NA,
             dx = NA,
             dy = NA)
    else
    {
        if (is.null(breaks)) {
            nint <- round(log2(length(x)) + 1)
            breaks <-
                if (equal.widths) do.breaks(range(x), nint)
                else quantile(x, 0:nint/nint)
        }
        h <- hist(x, breaks = breaks, plot = FALSE, ...)
        y <-
            if (type == "count") h$counts
            else if (type == "percent") 100 * h$counts/length(x)
            else h$intensities
        xlim <- range(x)
        lbreak <- max(xlim[1], breaks[breaks<=xlim[1]])
        ubreak <- min(xlim[2], breaks[breaks>=xlim[2]])

        list(xlim = range(x, lbreak, ubreak),
             ylim = range(0,y),
             dx = 1,
             dy = 1)
    }
}









panel.histogram <- function(x,
                            breaks,
                            equal.widths = TRUE,
                            type = "density",
                            col = bar.fill$col,
                            ...)
{
    grid.lines(x = c(0.05, 0.95),
               y = unit(c(0,0),"native"),
               default.units = "npc")
        
    if (length(x)>0) {
        bar.fill  <- trellis.par.get("bar.fill")

        if (is.null(breaks)) {

            nint <- round(log2(length(x)) + 1)
            breaks <-
                if (equal.widths) do.breaks(range(x), nint)
                else quantile(x, 0:nint/nint)

        }

        h <- hist(x, breaks = breaks, plot = FALSE, ...)
        y <-
            if (type == "count") h$counts
            else if (type == "percent") 100 * h$counts/length(x)
            else h$intensities

        nb <- length(breaks)
        if (nb != (length(y)+1)) warning("something is probably wrong")

        if (nb>1) {
            for(i in 1:(nb-1))
                if (y[i]>0) {
                    grid.rect(gp = gpar(fill = col),
                              x = breaks[i],
                              y = 0,
                              height = y[i],
                              width = breaks[i+1]-breaks[i],
                              just = c("left", "bottom"),
                              default.units = "native")
                }
        }
    }
}










histogram <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.histogram",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             type = c("percent", "count", "density"),
             nint = if (is.factor(x)) length(levels(x))
             else round(log2(length(x)) + 1),
             endpoints = extend.limits(range(x[!is.na(x)]), prop = 0.04),
             breaks = if (is.factor(x)) seq(0.5, length = length(levels(x))+1)
             else do.breaks(endpoints, nint),
             equal.widths = TRUE,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ## dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing

    formname <- deparse(substitute(formula))
    formula <- eval(substitute(formula), data, parent.frame())

    form <-
        if (inherits(formula, "formula"))
            latticeParseFormula(formula, data)
        else {
            if (!is.numeric(formula)) stop("invalid formula")
            else {
                list(left = NULL,
                     right = formula,
                     condition = NULL,
                     left.name = "",
                     right.name = formname)
            }
        }

    cond <- form$condition
    number.of.cond <- length(cond)
    x <- form$right
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    if(subscripts) subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    if (subscripts) subscr <- subscr[subset, drop = TRUE]
    
    if (missing(xlab)) xlab <- form$right.name
    if (missing(ylab)) ylab <- TRUE

    if(!(is.numeric(x) || is.factor(x)))
        warning("x should be numeric")
    x <- as.numeric(x)
    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))
                          

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- form$right.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo,
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:
    
    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)
        
        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        warning("Can't have log Y-scale")
        have.ylog <- FALSE
        foo$y.scales$log <- FALSE
    }

    if ((have.xlog || is.null(breaks) ||
        length(unique(round(diff(breaks)))) != 1) &&
        !missing(type))
        type <- "density"
    else type <- match.arg(type)

    if (is.logical(foo$ylab$label)) foo$ylab$label <- 
        if (type == "count") "Count"
        else if (type == "percent") "Percent of Total"
        else "Density"

        ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    ## equal.widths <- eval(equal.widths, data, parent.frame()) #keep this way ?
    foo$panel.args.common <- c(list(breaks = breaks,
                                    type = type,
                                    equal.widths = equal.widths), dots)
    if (subscripts) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout

    nplots <- plots.per.page * number.of.pages
    
    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    ##if (any(id)) {
                    foo$panel.args[[panel.number]] <-
                        list(x = x[id])
                    if (subscripts)
                        foo$panel.args[[panel.number]]$subscripts <-
                            subscr[id]
                    ##}
                    ##else
                    ##    foo$panel.args[[panel.number]] <-FALSE
                    
                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.histogram,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}



### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA




prepanel.default.levelplot <-
    function(x, y, wx, wy, subscripts, ...)
{
    xlim <- range(x[subscripts] + wx[subscripts]/2,
                  x[subscripts] - wx[subscripts]/2)
    ylim <- range(y[subscripts] + wy[subscripts]/2,
                  y[subscripts] - wy[subscripts]/2)
    list(xlim = extend.limits(xlim, prop = -0.0614),
         ylim = extend.limits(ylim, prop = -0.0614),
         dx = 1, dy = 1)
}




# index2seq <- function(x) {
#     print(x)
#     n <- length(x)
#     indices <- (1:n)[x==2]
#     ans <- rep(0, n)
#     for (i in indices)
#         ans[1:i] <- ans[1:i] + 1
#     ans
# }




panel.levelplot <-
    function(x, y, z, wx, wy, zcol, col.regions, subscripts,
             at = mean(z), labels = NULL, contour = TRUE, region = TRUE,
             col = add.line$col,
             lty = add.line$lty,
             lwd = add.line$lwd,
             ...)
{
    ## z not really needed here, but probably would be for contourplot
    if (any(subscripts)) {
        if (region) {
            for (i in seq(along = col.regions)) {
                ok <- (zcol[subscripts]==i)
                grid.rect(x = x[subscripts][ok],
                          y = y[subscripts][ok],
                          width = wx[subscripts][ok],
                          height = wy[subscripts][ok],
                          default.units = "native",
                          gp = gpar(fill=col.regions[i], col = NULL))
            }
        }
        if (contour) {
            add.line <- trellis.par.get("add.line")
            ux <- as.double(sort(unique(x[subscripts])))
            uy <- as.double(sort(unique(y[subscripts])))
            ord <- order(x[subscripts], y[subscripts])
            m <- z[subscripts][ord] + 10e-12
            for (i in seq(along = at)) {
                val <- .Call("cont", m, ux, uy, as.double(at[i]),
                             length(ux), length(uy), PACKAGE="lattice")
                if (is.null(labels))
                    lsegments(val[[1]], val[[2]], val[[3]], val[[4]],
                              col = col, lty = lty, lwd = lwd)
                else {
                    if (length(val[[1]]) > 3) {
                        ltext(lab = labels$lab[i],
                              x = .5 * (val[[1]][1]+val[[3]][1]),
                              y = .5 * (val[[2]][1]+val[[4]][1]))
                        lsegments(val[[1]][-(1:2)], val[[2]][-(1:2)],
                                  val[[3]][-(1:2)], val[[4]][-(1:2)])
                    }
                }
            }
        }
    }
}









contourplot <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.levelplot",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             cuts = 7,
             contour = TRUE,
             pretty = TRUE,
             region = FALSE,
             ...,
             subscripts = TRUE,
             subset = TRUE)

{
    ## m <- match.call(expand.dots = FALSE)
    dots <- list(...)
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    do.call("levelplot",
            c(list(formula = formula, data = data,
                   groups = groups, subset = subset,
                   panel = panel, prepanel = prepanel, strip = strip,
                   cuts = cuts, contour = contour,
                   pretty = pretty,
                   region = region),
              dots))
}







levelplot <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.levelplot",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             at,
             col.regions = trellis.par.get("regions")$col,
             colorkey = region,
             contour = FALSE,
             cuts = 15,
             labels = TRUE,
             pretty = FALSE,
             region = TRUE,
             ...,
             subscripts = TRUE,
             subset = TRUE)
{

    ##dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, z etc. and do some preprocessing

    formula <- eval(substitute(formula), data, parent.frame())
    form <-
        if (inherits(formula, "formula"))
            latticeParseFormula(formula, data, dim = 3)
        else {
            if (!is.matrix(formula)) stop("invalid formula")
            else {
                tmp <- expand.grid(1:nrow(formula), 1:ncol(formula))
                list(left = as.vector(formula),
                     right.x = tmp[[1]],
                     right.y = tmp[[2]],
                     condition = NULL,
                     left.name = "",
                     right.x.name = "", right.y.name = "")
            }
        }

    cond <- form$condition
    number.of.cond <- length(cond)
    z <- form$left
    x <- form$right.x
    y <- form$right.y

    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    y <- y[subset, drop = TRUE]
    z <- z[subset, drop = TRUE]
    subscr <- subscr[subset, drop = TRUE]

    if (missing(xlab)) xlab <- form$right.x.name
    if (missing(ylab)) ylab <- form$right.y.name

    if(!(is.numeric(x) && is.numeric(y) && is.numeric(z)))
        warning("x, y and z should be numeric")
    x <- as.numeric(x)
    y <- as.numeric(y)
    z <- as.numeric(z)

    zrng <- extend.limits(range(z[!is.na(z)]))
    if (missing(at))
        at <-
            if (pretty) lpretty(zrng, cuts)
            else seq(zrng[1], zrng[2], length = cuts+2)
    

    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))
                          

    ## Processing the labels argument
    if (is.logical(labels) && !labels) labels <- NULL
    else {
        if (is.logical(labels)) labels <- format(at)
        text <- trellis.par.get("par.xlab.text") # something better ?
        if (is.null(text)) text <- list(cex = 1, col = "black", font = 1, rot = 0)
        labels <- list(label = if (is.list(labels)) labels[[1]] else labels,
                       col = text$col, rot = text$rot,
                       cex = text$cex, font = text$font)
        if (is.list(labels)) labels[names(labels)] <- labels
        if (!is.character(labels$label))
            labels$label <- format(at)
    }

    
    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- form$right.name
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <- form$left.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo,
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        ylog <- foo$y.scales$log
        ybase <-
            if (is.logical(ylog)) 10
            else if (is.numeric(ylog)) ylog
            else if (ylog == "e") exp(1)

        y <- log(y, ybase)
        if (have.ylim) ylim <- log(ylim, ybase)
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))

    id.na <- is.na(x)|is.na(y)|is.na(z)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    ## Most levelplot/contourplot specific code here


    ## region
    numcol <- length(at)-1
    numcol.r <- length(col.regions)

    col.regions <-
        if (numcol.r <= numcol)
            rep(col.regions, length = numcol)
        else col.regions[floor(1+(1:numcol-1)*(numcol.r-1)/(numcol-1))]
    
    if (is.logical(colorkey)) {
        if (colorkey) foo$colorkey <-
            list(space = "right", col = col.regions,
                 at = at, tick.number = 7)
    }
    else if (is.list(colorkey)) {
        foo$colorkey <- colorkey
        if (is.null(foo$colorkey$col)) foo$colorkey$col <- col.regions
        if (is.null(foo$colorkey$at)) foo$colorkey$at <- at
    }

    
    ## I'm going to create vectors parallel to x y etc which would
    ## give the widths and heights of the rectangles for each point.
    ## My algo works only when the x's and y's are really evaluated
    ## on a grid, that is, there is no numerical error. Splus also
    ## doesn't work (in any meningful way, at least) in such cases,
    ## but behaviour would be dissimilar in that case.

    ux <- sort(unique(x[!is.na(x)]))
    dux <- diff(ux)
    wux <- .5 * (c(dux[1], dux) + c(dux, dux[length(dux)]))
    ##wx <- wux[match(x[!is.na(x)], ux)]
    wx <- wux[match(x, ux)]
    uy <- sort(unique(y[!is.na(y)]))
    duy <- diff(uy)
    wuy <- .5 * (c(duy[1], duy) + c(duy, duy[length(duy)]))
    ##wy <- wuy[match(y[!is.na(y)], uy)]
    wy <- wuy[match(y, uy)]

    zcol <- numeric(length(z))
    for (i in seq(along=col.regions))
        zcol[!id.na & z>=at[i] & z<at[i+1]] <- i

    foo$panel.args.common <-
        c(list(x=x, y=y, z=z, at=at, labels=labels,
               region = region, contour = contour,
               wx=wx, wy=wy, zcol=zcol, col.regions=col.regions),
          dots)

    if (!is.null(groups)) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1, number.of.cond)
    panel.number <- 1 # this is a counter for panel number

    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    foo$panel.args[[panel.number]] <- 
                        list(subscripts = subscr[id])

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.levelplot,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}











### Copyright 2001-2002 Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA


lpretty <- function(x, ...) { 
    eps <- 1e-10
    at <- pretty(x, ...)
    at <- ifelse(abs(at-round(at, 3))<eps, round(at, 3), at)
}


oneway <-
    function(formula, data, location = mean,
             spread = function(x) sqrt(var(x)))
{
    if(missing(data)) data <- sys.frame(sys.parent())
    form <- latticeParseFormula(formula, data)
    y <- form$left
    x <- form$right
    if (!is.shingle(x)) x <- as.factor(x)
    is.f.x <- is.factor(x)
    num.l.x <- nlevels(x) 
    foo <- list()
    if (is.f.x) {
        foo$location <-
            if (is.function(location)) as.vector(tapply(X=y, INDEX=list(x), FUN = location))
            else rep(location, num.l.x)
        foo$spread <- 
            if (is.function(spread)) as.vector(tapply(X=y, INDEX=list(x), FUN = spread))
            else rep(spread, num.l.x)
        foo$fitted.values <- numeric(length(y))
        sc <- numeric(length(y))
        for (i in seq(along = y)){
            foo$fitted.values[i] <- foo$location[as.numeric(x)[i]]
            sc[i] <- foo$spread[as.numeric(x)[i]]
        }
        foo$residuals <- y - foo$fitted.values
        foo$scaled.residuals <- foo$residuals/sc
    }
    else stop("x must be (coercible to be) a factor")
    foo
}


do.breaks  <- function(endpoints, nint)
{
    if (length(endpoints)!=2) stop("error")
    endpoints[1] + diff(endpoints) * 0:nint / nint
}


## This converts character to factor, numeric to shingle, and
## in addition, takes subsets
as.factorOrShingle <- function(x, subset = TRUE, drop = FALSE)
{
    if (is.character(x)) as.factor(x)[subset, drop = drop]
    else if (is.numeric(x)) as.shingle(x)[subset, drop = drop]
    else x[subset, drop = drop]
}



"[.shingle" <-
    function(x, subset, drop = FALSE)
{
    if (!is.shingle(x)) stop("x must be a shingle")
    ans <- as.numeric(x)[subset]
    attr(ans, "levels") <- levels(x)
    class(attr(ans, "levels")) <- "shingleLevel"
    if (drop) {
        xlvs <- levels(ans)
        dl <- logical(nlevels(ans))
        for (i in seq(along=dl))
            dl[i] <- any( ans >= xlvs[[i]][1] & ans <= xlvs[[i]][2] )
        attr(ans, "levels") <- xlvs[dl]
        class(attr(ans, "levels")) <- "shingleLevel"
    }
    class(ans) <- "shingle"
    ans
}



Rows <- function(x, which)
{
    for (i in seq(along = x)) x[[i]] <- x[[i]][which]
    x
}
## S-Plus trellis function needed for nlme.



make.list.from.intervals <- function(x)
{
    if (ncol(x)!=2) stop("x must be matrix with 2 columns")
    if (nrow(x)<1) stop("x must be matrix with at least 1 row")
    ans <- as.list(1:nrow(x))
    for (i in 1:nrow(x))
        ans[[i]] <- x[i,]
    ans
}



equal.count <-
  function(x, ...)
{
    attr(x, "levels") <- make.list.from.intervals(co.intervals(x,...))
    class(attr(x, "levels")) <- "shingleLevel"
    class(x) <- "shingle"
    x
}



shingle <-
    function(x, intervals=sort(unique(x)))
{
    if (ncol(as.matrix(intervals))==1)
        intervals <- cbind(intervals, intervals)
    else if (ncol(as.matrix(intervals)) > 2)
        stop("bad value of 'intervals'")
    attr(x, "levels") <- make.list.from.intervals(intervals)
    class(attr(x, "levels")) <- "shingleLevel"
    class(x) <- "shingle"
    x
}


as.data.frame.shingle <- as.data.frame.factor

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


as.shingle <-
    function(x) if (is.shingle(x)) x else shingle(x)



summary.shingle <- function(object, ...) print.shingle(object, ...)


print.shingleLevel <-
    function(x, ...) {
        print(do.call("rbind", x))
        invisible(x)
    }

print.shingle <- function(x, ...) {
    cat("\nData:\n")
    print(as.numeric(x))
    l <- levels(x)
    n <- nlevels(x)
    if (n<1) cat("\nno intervals\n")
    else {
        int <- data.frame(min = numeric(n), max = numeric(n), count = numeric(n))
        for (i in 1:n) {
            int$min[i] <- l[[i]][1]
            int$max[i] <- l[[i]][2]
            int$count[i] <- length(x[x>=l[[i]][1] & x<=l[[i]][2]])
        }
        cat("\nIntervals:\n")
        print(int)
        olap <- numeric(n-1)
        if (n>2)
            for (i in 1:(n-1))
                olap[i] <- length(x[ x>=l[[i]][1] & x<=l[[i]][2] &
                                    x>=l[[i+1]][1] & x<=l[[i+1]][2]])
        cat("\nOvrlap between adjacent intervals:\n")
        print(olap)
    }
    invisible(x)
}





strip.default <-
    function(which.given,
             which.panel,
             var.name,
             factor.levels,
             shingle.intervals,
             strip.names = c(FALSE, TRUE),
             style = 1,
             bg = trellis.par.get("strip.background")$col[which.given],
             fg = trellis.par.get("strip.shingle")$col[which.given],
             par.strip.text = trellis.par.get("add.text"))
{
    name <- var.name[which.given]
    level <- which.panel[which.given]
    strip.names <- rep(strip.names, length = 2)
    
    if (is.null(factor.levels)) { # means this is a  shingle, as opposed to a factor
        if (is.null(shingle.intervals)) stop("both factor.levels and shingle.intervals cannot be NULL")
        strip.names <- strip.names[2]
        grid.rect(gp = gpar(fill=bg))
        t <- range(shingle.intervals)
        r <- (range(shingle.intervals[level,])-t[1])/diff(t)
        grid.rect(x = unit(r%*%c(.5,.5),"npc"), width = unit(diff(r),"npc"),
                  gp = gpar(col=fg, fill=fg))
        if (strip.names) grid.text(label = name,
                                   gp = gpar(col = par.strip.text$col,
                                   font = par.strip.text$font,
                                   fontsize = par.strip.text$cex *
                                   current.viewport()$gp$fontsize))
        
        grid.rect()
    }
    else if (is.null(shingle.intervals)) { # factor
        strip.names <- strip.names[1]
        x <- factor.levels
        num <- length(x)
        if (style == 1) {
            grid.rect(gp = gpar(fill=bg))
            grid.text(label = paste(if(strip.names)
                      paste(name,": ") else "", x[level], sep = ""),
                      gp = gpar(col = par.strip.text$col,
                      font = par.strip.text$font,
                      fontsize = par.strip.text$cex *
                      current.viewport()$gp$fontsize))
            grid.rect()
        }
        else if (style == 2) {
            grid.rect(x = unit((2*level-1)/(2*num), "npc"),
                      width = unit(1/num, "npc"),
                      gp = gpar(fill=fg, col = NULL))
            grid.text(label=x,
                      x = (2*1:num-1)/(2*num),
                      gp = gpar(col = par.strip.text$col,
                      font = par.strip.text$font,
                      fontsize = par.strip.text$cex *
                      current.viewport()$gp$fontsize))
            grid.rect()
        }
        else if (style == 3){
            grid.rect(gp = gpar(fill=bg))
            grid.rect(x = unit((2*level-1)/(2*num), "npc"),
                      width = unit(1/num, "npc"),
                      gp = gpar(fill=fg, col = NULL))
            grid.text(label = paste(if(strip.names)
                      paste(name,": ") else "", x[level], sep = ""),
                      gp = gpar(col = par.strip.text$col, 
                      font = par.strip.text$font,
                      fontsize = par.strip.text$cex *
                      current.viewport()$gp$fontsize))
            grid.rect()
        }
        else if(style == 4){
            grid.rect(gp = gpar(fill=bg))
            grid.rect(x = unit((2*level-1)/(2*num), "npc"),
                      width = unit(1/num, "npc"),
                      gp = gpar(col=NULL, fill=fg))
            grid.text(label=x,
                      x = (2* 1:num - 1)/(2*num),   #using default.units
                      gp = gpar(col = par.strip.text$col, 
                      font = par.strip.text$font,
                      fontsize = par.strip.text$cex *
                      current.viewport()$gp$fontsize))
            grid.rect()
        }
        else if(style >= 5){
            grid.rect(gp = gpar(fill=bg))
            grid.text(label=x[level],
                      x = (2* level - 1)/(2*num),   #using default.units
                      gp = gpar(col = par.strip.text$col, 
                      font = par.strip.text$font,
                      fontsize = par.strip.text$cex *
                      current.viewport()$gp$fontsize))
            grid.rect()
        }
    }
}







lsegments <-
    function(x0, y0, x1, y1, 
             col = add.line$col,
             lty = add.line$lty,
             lwd = add.line$lwd, ...)
{
    add.line <- trellis.par.get("add.line")
    ml <- max(length(x0), length(x1), length(y0), length(y1))
    x0 <- rep(x0, length = ml)
    x1 <- rep(x1, length = ml)
    y0 <- rep(y0, length = ml)
    y1 <- rep(y1, length = ml)

    grid.segments(x0 = x0, x1 = x1,
                  y0 = y0, y1 = y1,
                  gp = gpar(lty=lty,
                  col=col, lwd=lwd),
                  default.units="native")
}


larrows <-
    function(x0, y0, x1, y1, angle = 30, code = 2, length = NULL, proportion = .05, ...) 
{
    if (!is.null(length)) warning("length not implemented in larrows, use proportion instead")
    ##warning("larrows not corectly implemented yet")
    angle <- angle / 180 * pi
    start <- rbind(x0, y0)
    end <- rbind(x1, y1)
    v.forward <- end - start
    v.backward <- start - end
    lsegments(x0, y0, x1, y1, ...)
    
    if (code %in% c(1,3)) { # arrow at starting point
        edge.1 <- proportion * 
            matrix( c(cos(angle), -sin(angle), sin(angle), cos(angle)), 2, 2) %*% v.forward
        edge.2 <- proportion *
            matrix( c(cos(-angle), -sin(-angle), sin(-angle), cos(-angle)), 2, 2) %*% v.forward
        lsegments(x0, y0, x0 + edge.1[1,], y0 + edge.1[2,], ...)
        lsegments(x0, y0, x0 + edge.2[1,], y0 + edge.2[2,], ...)
    }
    if (code %in% c(2,3)) { # arrow at ending point
        edge.1 <- proportion * 
            matrix( c(cos(angle), -sin(angle), sin(angle), cos(angle)), 2, 2) %*% v.backward
        edge.2 <- proportion *
            matrix( c(cos(-angle), -sin(-angle), sin(-angle), cos(-angle)), 2, 2) %*% v.backward
        lsegments(x1, y1, x1 + edge.1[1,], y1 + edge.1[2,], ...)
        lsegments(x1, y1, x1 + edge.2[1,], y1 + edge.2[2,], ...)
    }
}



ltext <-
    function(x, y = NULL, labels = seq(along = x),
             col = add.text$col,
             cex = add.text$cex,
             srt = 0,
             font = 1,
             adj = .5, ...)
{
    add.text <- trellis.par.get("add.text")
    xy <- xy.coords(x, y)
    grid.text(label = as.character(labels), x = xy$x, y = xy$y,
              gp = gpar(col = col, font = font,
              fontsize = cex * 10),
              just = c(if (adj == 0) "left"
              else if (adj == 1) "right" else "centre", "centre"),
              rot = srt,
              default.units = "native")
}





llines <-
    function(x, y = NULL, type = "l",
             col = plot.line$col,
             lty = plot.line$lty,
             lwd = plot.line$lwd, ...)
{
    plot.line <- trellis.par.get("plot.line")
    lplot.xy(xy.coords(x, y), type = type,
             col = col, lty = lty, lwd = lwd, ...)
}




lpoints <-
    function(x, y = NULL, type = "p",
             col = plot.symbol$col,
             pch = plot.symbol$pch,
             cex = plot.symbol$cex, ...)
{
    plot.symbol <- trellis.par.get("plot.symbol")
    lplot.xy(xy.coords(x, y), type = type,
             col = col, pch = pch, cex = cex, ...)
}






lplot.xy <-
    function(xy, type, pch = 1, lty = 1, col = 1, cex = 1, lwd = 1, font = 1, ...)
{
    x <- xy$x
    y <- xy$y

    if (type %in% c("l", "o", "b", "c"))
        grid.lines(x=x, y=y, gp = gpar(lty=lty, col=col, lwd=lwd),
                   default.units="native")
    
    if (type %in% c("p", "o", "b", "c"))
        if (is.character(pch))
            grid.text(lab = rep(pch, length = length(x)),
                      x = x, y = y,
                      gp = gpar(col = col, fontsize = cex * 10),
                      default.units="native")
        else
            grid.points(x = x, y = y, size = unit(cex * 2.5, "mm"),
                        gp = gpar(col = col),
                        pch = pch, 
                        default.units="native")

    if (type %in% c("s", "S")) {
        ord <- order(x)
        n <- length(x)
        xx <- numeric(2*n-1)
        yy <- numeric(2*n-1)

        xx[2*1:n-1] <- x[ord]
        yy[2*1:n-1] <- y[ord]
        xx[2*1:(n-1)] <- x[ord][if (type=="s") -1 else -n]
        yy[2*1:(n-1)] <- y[ord][if (type=="s") -n else -1]
        grid.lines(x=xx, y=yy,
                   gp = gpar(lty=lty, col=col, lwd=lwd),
                   default.units="native")
    }

    if (type == "h") {

        ylim <- current.viewport()$yscale

        zero <-
            if (ylim[1] > 0) ylim[1]
            else if (ylim[2] < 0) ylim[2]
            else 0

        ##print(zero) ?
        ##print(x) 
        for (i in seq(along=x))
            grid.lines(x=rep(x[i],2), y=c(y[i], zero),
                       gp = gpar(lty=lty, col=col, lwd=lwd),
                       default.units="native")
    }
}









### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA

## the foll functions don't do much error checking yet



panel.abline <-
    function(a, b = NULL, h = numeric(0), v = numeric(0),
             col = add.line$col, lty = add.line$lty,
             lwd = add.line$lwd, ...)
{
    add.line <- trellis.par.get("add.line")
    
    if (!missing(a)) {
        if (inherits(a,"lm")) {
            coeff <- coef(a)
        }
        else if (!is.null(coef(a))) coeff <- coef(a)  # ????
        else coeff <- c(a,b)

        if (length(coeff)==1) coeff <- c(0, coeff)
        
        if (coeff[2]==0) h <- c(h, coeff[1])
        else {
            xx <- current.viewport()$xscale
            yy <- current.viewport()$yscale
            
            x <- numeric(0)
            y <- numeric(0)
            ll <- function(i, j, k, l)
                (yy[j]-coeff[1]-coeff[2]*xx[i]) *
                    (yy[l]-coeff[1]-coeff[2]*xx[k])
            
            if (ll(1,1,2,1)<=0) {
                y <- c(y, yy[1])
                x <- c(x, (yy[1]-coeff[1])/coeff[2])
            }
            
            if (ll(2,1,2,2)<=0) {
                x <- c(x, xx[2])
                y <- c(y, coeff[1] + coeff[2] * xx[2])
            }
            
            if (ll(2,2,1,2)<=0) {
                y <- c(y, yy[2])
                x <- c(x, (yy[2]-coeff[1])/coeff[2])
            }
            
            if (ll(1,2,1,1)<=0) {
                x <- c(x, xx[1])
                y <- c(y, coeff[1] + coeff[2] * xx[1])
            }
            
            if (length(x)>0)
                grid.lines(x=x, y = y, default.units="native",
                           gp = gpar(col=col, lty=lty, lwd=lwd))
        }
    }
    
    
    for(i in seq(along=h))
        grid.lines(y=rep(h[i],2), default.units="native", gp = gpar(col=col,lty=lty,lwd=lwd))

    for(i in seq(along=v))
        grid.lines(x=rep(v[i],2), default.units="native", gp = gpar(col=col,lty=lty,lwd=lwd))
    
}








panel.fill <-
    function(col = trellis.par.get("background")$col, ...)
{
    grid.rect(gp=gpar(fill=col))
}












panel.grid <-
    function(h=3, v=3, col=reference.line$col, lty=reference.line$lty,
             lwd=reference.line$lwd, ...)
{
    reference.line <- trellis.par.get("reference.line")

    if (h>0)
        for(i in 1:h)
            grid.lines(y=rep(i/(h+1),2),
                       gp = gpar(col = col, lty = lty, lwd = lwd),
                       default.units="npc")

    if (v>0)
        for(i in 1:v)
            grid.lines(x=rep(i/(v+1),2),
                       gp = gpar(col = col, lty = lty, lwd = lwd),
                       default.units="npc")


    ## Cheating here a bit for h=-1, v=-1. Can't think of any neat way to
    ## get the actual `at' values of the panel (Can pass it in though)

    if (h<0)
    {
        scale <- current.viewport()$yscale
        at <- lpretty(scale)
        at <- at[at>scale[1] & at < scale[2]]
        for(i in seq(along=at))
            grid.lines(y=rep(at[i],2), default.units="native",
                       gp = gpar(col = col, lty = lty, lwd = lwd))
    }
    if (v<0)
    {
        scale <- current.viewport()$xscale
        at <- lpretty(scale)
        at <- at[at>scale[1] & at < scale[2]]
        for(i in seq(along=at))
            grid.lines(x=rep(at[i],2), default.units="native",
                       gp = gpar(col = col, lty = lty, lwd = lwd))
    }
}





panel.lmline <-
    function(x, y, ...) if (length(x)>0) panel.abline(lm(y ~ x), ...) 



prepanel.lmline <-
    function(x, y, ...)
{
    if (length(x)>0) {
        coeff <- coef(lm(y~x))
        tem <- coeff[1] + coeff[2] * range(x)
        list(xlim=range(x), ylim=range(y,tem), 
             dx=diff(range(x)), dy=diff(tem))         
    }
    else list(xlim=c(NA,NA), ylim=c(NA,NA), dx=NA, dy=NA)
}


panel.loess <-
    function(x, y, span = 2/3, degree = 1,
             family = c("symmetric", "gaussian"),
             evaluation = 50,
             lwd = add.line$lwd, lty = add.line$lty,
             col = add.line$col, ...)
{
    if (length(x)>0) {
        add.line <- trellis.par.get("add.line")
        
        smooth <- loess.smooth(x, y, span = span, family = family,
                               degree = degree, evaluation = evaluation)
        grid.lines(x=smooth$x, y=smooth$y, default.units = "native",
                   gp = gpar(col = col, lty = lty, lwd = lwd))
    }
}


prepanel.loess <-
    function(x, y, span = 2/3, degree = 1,
             family = c("symmetric", "gaussian"),
             evaluation = 50,
             lwd = add.line$lwd, lty = add.line$lty,
             col = add.line$col, ...)
{
    if (length(x)>0) {
        add.line <- trellis.par.get("add.line")
        
        smooth <- loess.smooth(x, y, span = span, family = family,
                               degree = degree, evaluation = evaluation)
        list(xlim = range(x,smooth$x),
             ylim = range(y,smooth$y),
             dx = diff(smooth$x),
             dy = diff(smooth$y))
    }
    else list(xlim=c(0,1), ylim=c(0,1), dx=1, dy=1)
}



# panel.smooth <-
#     function(x, y, span = 2/3, degree = 1, zero.line = FALSE,
#              family = c("symmetric", "gaussian"),
#              evaluation = 50,
#              lwd = add.line$lwd, lty = add.line$lty,
#              col = add.line$col, ...)
# {
#     if (zero.line) abline(h=0, ...)
#     panel.loess(x, y, span = span, family = family,
#                 degree = degree, evaluation = evaluation, ...)
#     panel.xyplot(x, ,y, ...)
# }
## base R function exists




panel.superpose <-
    function(x, y = NULL, subscripts, groups,
             panel.groups = "panel.xyplot",
             col,
             col.line = superpose.line$col,
             col.symbol = superpose.symbol$col,
             pch = superpose.symbol$pch,
             cex = superpose.symbol$cex, 
             lty = superpose.line$lty,
             lwd = superpose.line$lwd,
             ...)
{
    if (length(x)>0) {

        if (!missing(col)) {
            if (missing(col.line)) col.line <- col
            if (missing(col.symbol)) col.symbol <- col
        }

        superpose.symbol <- trellis.par.get("superpose.symbol")
        superpose.line <- trellis.par.get("superpose.line")

        x <- as.numeric(x)
        if (!is.null(y)) y <- as.numeric(y)

        vals <- sort(unique(groups))
        nvals <- length(vals)
        col.line <- rep(col.line, length=nvals)
        col.symbol <- rep(col.symbol, length=nvals)
        pch <- rep(pch, length=nvals)
        lty <- rep(lty, length=nvals)
        lwd <- rep(lwd, length=nvals)
        cex <- rep(cex, length=nvals)

        panel.groups <- 
            if (is.function(panel.groups)) panel.groups
            else if (is.character(panel.groups)) get(panel.groups)
            else eval(panel.groups)

        for (i in seq(along=vals)) {
            id <- (groups[subscripts] == vals[i])
            if (any(id)) {
                args <- list(x=x[id], 
                             pch = pch[i], cex = cex[i],
                             col.line = col.line[i],
                             col.symbol = col.symbol[i],
                             lty = lty[i],
                             lwd = lwd[i], ...)
                if (!is.null(y)) args$y=y[id]

                do.call("panel.groups", args)
            }
        }
    }
}












panel.superpose.2 <- 
    function (x, y, subscripts, groups, col, col.line = superpose.line$col,
              col.symbol = superpose.symbol$col, pch = superpose.symbol$pch,
              cex = superpose.symbol$cex, lty = superpose.line$lty,
              lwd = superpose.line$lwd, type="p", ...)
{
    
    ##   `panel.superpose.2' :  This is a version of the 'panel.superpose'
    ##   Trellis panel function that allows the plot `type' to change between
    ##   superimposed (overlayed) data sets.  See the `panel.xyplot' function
    ##   for details on the `type' option which is usually a single character,
    ##   but here is a character vector with each element specifying the
    ##   plot style of each subsequently-overlayed plot.
    ##                        ---  Neil Klepeis, 26-Dec-2001
    
    if (length(x) > 0) {
        if (!missing(col)) {
            if (missing(col.line))
                col.line <- col
            if (missing(col.symbol))
                col.symbol <- col
        }
        superpose.symbol <- trellis.par.get("superpose.symbol")
        superpose.line <- trellis.par.get("superpose.line")
        x <- as.numeric(x)
        y <- as.numeric(y)
        vals <- sort(unique(groups))
        nvals <- length(vals)
        col.line <- rep(col.line, length = nvals)
        col.symbol <- rep(col.symbol, length = nvals)
        pch <- rep(pch, length = nvals)
        lty <- rep(lty, length = nvals)
        lwd <- rep(lwd, length = nvals)
        cex <- rep(cex, length = nvals)
        type <- rep(type, length = nvals)      # new line here
        for (i in seq(along = vals)) {
            id <- (groups[subscripts] == vals[i])
            if (any(id))
                panel.xyplot(x = x[id], y = y[id], pch = pch[i],
                  cex = cex[i], col.line = col.line[i], col.symbol = col.symbol[i],
                  lty = lty[i], lwd = lwd[i], type=type[i], ...)
        }
    }
}







panel.linejoin <-
    function(x, y, fun = mean,
             horizontal = TRUE,
             lwd = reference.line$lwd,
             lty = reference.line$lty,
             col = reference.line$col,
             col.line,
             ...)
{
    reference.line = trellis.par.get("reference.line")
    if (missing(col.line)) col.line <- col
    if (horizontal) {
        vals <- unique(sort(y))
        yy <- seq(along = vals)
        xx <- numeric(length(yy))
        for (i in yy)
            xx[i] <- fun(x[y == vals[i]])
        llines(xx, yy, col = col.line, lty = lty, lwd = lwd, ...)
    }
    else {
        vals <- unique(sort(x))
        xx <- seq(along = vals)
        yy <- numeric(length(xx))
        for (i in xx)
            yy[i] <- fun(y[x == vals[i]])
        llines(xx, yy, col = col.line, lty = lty, lwd = lwd, ...)
     }
}



panel.mathdensity <-
    function(dmath = dnorm,
             args = list(mean = 0, sd = 1),
             n = 50,
             col = reference.line$col,
             lwd = reference.line$lwd, ...)
{

    reference.line <- trellis.par.get("reference.line")
    x <- do.breaks(endpoints = current.viewport()$xscale,
                   nint = n)
    y <- do.call("dmath", c(list(x = x),args))
    panel.xyplot(x = x, y = y, type = "l", col = col, lwd = lwd, ...)
    
}


### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



prepanel.default.parallel <-
    function(x, y, type, ...)
{
    list(xlim = c(0,1),
         ylim = c(0,1),
         dx = 1,
         dy = 1)
}



panel.parallel <- function(z, subscripts,
                           col=superpose.line$col,
                           lwd=superpose.line$lwd,
                           lty=superpose.line$lty, ...)
{
    superpose.line <- trellis.par.get("superpose.line")
    reference.line <- trellis.par.get("reference.line")

    n.r <- ncol(z)
    n.c <- length(subscripts)
    col <- rep(col, length=n.c)
    lty <- rep(lty, length=n.c)
    lwd <- rep(lwd, length=n.c)

    llim <- numeric(n.r)
    ulim <- numeric(n.r)
    dif <- numeric(n.r)
    if(n.r>0)
        for(i in 1:n.r) {
            grid.lines(x = c(0,1), y = c(i,i),
                       default.units = "native",
                       gp = gpar(col = reference.line$col,
                       lwd = reference.line$lwd,
                       lty = reference.line$lty))
            llim[i] <- range(z[,i])[1]
            ulim[i] <- range(z[,i])[2]
            dif[i] <- ulim[i] - llim[i]
        }
   

    for (i in seq(along=subscripts))
    {
        x <- (as.numeric(z[subscripts[i],,])-llim)/dif
        grid.lines(x = x,
                   y=1:n.r, 
                   gp = gpar(col=col[i], lty=lty[i], lwd=lwd[i]),
                   default.units="native")
    }
    
}






parallel <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             between = list(x = 0.5, y = 0.5),
             layout = NULL,
             panel = "panel.parallel",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab = NULL,
             xlim,
             ylab = NULL,
             ylim,
             varnames,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ## dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing
    
    form <- latticeParseFormula(formula, data)
    cond <- form$condition
    number.of.cond <- length(cond)
    x <- as.data.frame(form$right)
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, nrow(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }
    if (!missing(varnames)) colnames(x) <-
        eval(substitute(varnames), data, parent.frame())

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    subscr <- seq(along=x[,1])
    x <- x[subset,, drop = TRUE]
    subscr <- subscr[subset, drop = TRUE]
    
    ##if(!(is.numeric(x) && is.numeric(y)))
    ##    warning("Both x and y should be numeric")    WHAT ?


    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          between = between,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- "Scatter Plot Matrix"
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab <- NULL

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## overriding at and labels, maybe not necessary
    
    ## scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    if (is.null(scales$alternating)) {
        if (is.null(scales$y)) scales$y <- list(alternating = FALSE)
        else if (is.null(scales$y$alternating)) scales$y$alternating <- FALSE
        ## bug if y="free" but who cares
    }
    foo <- c(foo, 
             do.call("construct.scales", scales))
    foo$x.scales$at <- c(0,1)
    foo$x.scales$labels <- c("Min","Max")
    foo$y.scales$at <- 1:ncol(x)
    foo$y.scales$labels <- colnames(x)
    

    ## Step 3: Decide if limits were specified in call:

    if (missing(xlim)) xlim <- extend.limits(c(0,1))
    if (missing(ylim)) ylim <- extend.limits(c(1,ncol(x)), prop = 0.03) 
    have.xlim <- TRUE
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- TRUE
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }
    
    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        foo$x.scales$log <- FALSE
        ## This is because No further changes will be
        ## necessary while printing since x-axes are not
        ## marked (many x axes)
    }
    if (have.ylog) {
        warning("cannot have log y-scale")
        foo$y.scales$log <- FALSE
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- FALSE
    for (j in 1:ncol(x))
        id.na <- id.na | is.na(x[,j])
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args


    foo$panel.args.common <-
        c(list(z = x, groups = groups), dots)

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    foo$panel.args[[panel.number]] <-
                        list(subscripts = subscr[id])

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.parallel,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}















### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



## not quite what it should be
plot.shingle <-
    function(x, col = bar.fill$col, aspect = "fill", ...)
{

    bar.fill <- trellis.par.get("bar.fill")
    foo <- list(call = match.call(),
                aspect.fill = aspect == "fill",
                aspect.ratio = if (is.numeric(aspect)) aspect else 1,
                as.table = FALSE,
                condlevels = "1",
                key = NULL,
                layout=c(1,1,1),
                page = NULL,
                panel = function(x, col) {
                    ## x is the list of intervals
                    num.l.y <- length(x)
                    if (num.l.y>0)
                        for(i in 1:num.l.y)
                            grid.rect(x = x[[i]] %*% c(.5,.5),
                                      y = i,
                                      width = diff(x[[i]]),
                                      height = .5,
                                      default.units = "native",
                                      gp = gpar(fill=col)) 
                },
                panel.args = list(list()),
                panel.args.common = list(x=levels(x), col = col),
                par.strip.text = trellis.par.get("add.text"),
                skip = FALSE,
                strip = FALSE,
                main = NULL,
                sub = NULL,
                xlab = list(label = "Range", col = "black", cex = 1, font =1),
                ylab = list(label = "Panel", col = "black", cex = 1, font =1),
                x.scales = 1,
                y.scales = 1,
                x.between = 0,
                y.between = 0,
                x.alternating = 1,
                y.alternating = 1,
                fontsize.normal = 10,
                fontsize.small = 8)
    
    num.l.y <- nlevels(x)
    foo$x.limits <- extend.limits(range(x, levels(x)))
    foo$y.limits <- extend.limits(c(1,num.l.y),
                                  length = .5+num.l.y)


    foo$x.scales <- list(relation = "same",
                         draw = TRUE,
                         alternating = 1,
                         at = FALSE,
                         labels = FALSE,
                         tck = 1,
                         font = 1,
                         col = FALSE,
                         log = FALSE,
                         cex = 1,
                         rot = FALSE,
                         tick.number = 5)
    
    foo$y.scales <- list(relation = "same",
                         draw = TRUE,
                         alternating = 1,
                         at = 1:num.l.y,
                         labels = FALSE,
                         tck = 1,
                         font = 1,
                         col = FALSE,
                         log = FALSE,
                         cex = 1,
                         rot = FALSE,
                         tick.number = num.l.y)
    
    class(foo) <- "trellis"
    foo
    
}






### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



draw.key <- function(key, draw = FALSE, vp = NULL)
{
    if (!is.list(key)) stop("key must be a list")
    
    max.length <- 0
    ## maximum of the `row-lengths' of the above
    ## components. There is some scope for confusion
    ## here, e.g., if col is specified in key as a
    ## length 6 vector, and then lines=list(lty=1:3),
    ## what should be the length of that lines column ?
    ## If 3, what happens if lines=list() ?
    ## (Strangely enough, S+ accepts lines=list()
    ## if col (etc) is NOT specified outside, but not
    ## if it is)
    
    process.key <-
        function(between = 2,
                 align = TRUE,
                 title = "",
                 background = trellis.par.get("background")$col,
                 border = NULL,
                 transparent = FALSE, 
                 columns = 1,
                 divide = 3,
                 between.columns = 3,
                 cex = 1,
                 cex.title = 1.5 * max(cex),
                 col = "black", 
                 lty = 1,
                 lwd = 1,
                 font = 1, 
                 pch = 8,
                 adj = 0,
                 type = "l", 
                 size = 5, 
                 angle = 0, 
                 density = -1,
                 ...)
        {
            list(between = between,
                 align = align,
                 title = title,
                 background = background,
                 border = border,
                 transparent = transparent, 
                 columns = columns,
                 divide = divide,
                 between.columns = between.columns,
                 cex = cex,
                 cex.title = cex.title,
                 col = col,
                 lty = lty,
                 lwd = lwd,
                 font = font, 
                 pch = pch,
                 adj = adj,
                 type = type, 
                 size = size, 
                 angle = angle, 
                 density = density,
                 ...)
        }

    key <- do.call("process.key", key)

    key.length <- length(key)
    key.names <- names(key)    # Need to update
    if (is.logical(key$border)) 
        key$border <-
            if (key$border) "black"
            else NULL

    components <- list()

    for(i in 1:key.length) {

        curname <- pmatch(key.names[i], c("text", "rectangles", "lines", "points"))

        if (is.na(curname)) {
            ;## do nothing
        }
        else if (curname == 1) { # "text"
            if (!(is.character(key[[i]][[1]])))
                stop("first component of text has to be vector of labels")
            pars <- list(labels = as.character(key[[i]][[1]]),
                         col = key$col,
                         cex = key$cex,
                         font = key$font)
            key[[i]][[1]] <- NULL
            pars[names(key[[i]])] <- key[[i]]

            tmplen <- length(pars$labels)
            for (j in 1:length(pars))
                pars[[j]] <- rep(pars[[j]], length = tmplen)

            max.length <- max(max.length, tmplen)
            components[[length(components)+1]] <-
                list(type = "text", pars = pars, length = tmplen)

        }
        else if (curname == 2) { # "rectangles"

            pars <- list(col = key$col,
                         size = key$size,
                         angle = key$angle,
                         density = key$density)
            
            pars[names(key[[i]])] <- key[[i]]

            tmplen <- max(unlist(lapply(pars,length)))
            max.length <- max(max.length, tmplen)
            components[[length(components)+1]] <-
                list(type = "rectangles", pars = pars)
            
        }
        else if (curname == 3) { # "lines"

            pars <- list(col = key$col,
                         size = key$size,
                         lty = key$lty,
                         lwd = key$lwd,
                         type = key$type)
                         
            pars[names(key[[i]])] <- key[[i]]

            tmplen <- max(unlist(lapply(pars,length)))
            max.length <- max(max.length, tmplen)
            components[[length(components)+1]] <-
                list(type = "lines", pars = pars)
            
        }
        else if (curname == 4) { # "points"

            pars <- list(col = key$col,
                         cex = key$cex,
                         pch = key$pch)
                         
            pars[names(key[[i]])] <- key[[i]]

            tmplen <- max(unlist(lapply(pars,length)))
            max.length <- max(max.length, tmplen)
            components[[length(components)+1]] <-
                list(type = "points", pars = pars)

        }
    }

    number.of.components <- length(components)
    ## number of components named one of "text",
    ## "lines", "rectangles" or "points"

    for (i in 1:number.of.components)
        if (components[[i]]$type != "text") {
            components[[i]]$pars <-
                lapply(components[[i]]$pars, rep, length = max.length)
            components[[i]]$length <- max.length
        }

    column.blocks <- key$columns
    rows.per.block <- ceiling(max.length/column.blocks)

    if (column.blocks > max.length) warning("not enough rows for columns")
    
    key$between <- rep(key$between, length = number.of.components)

    
    if (key$align) {

        ## Setting up the layout

        n.row <- rows.per.block + 1
        n.col <- column.blocks * (1 + 3 * number.of.components) - 1

        heights.x <- rep(1, n.row)
        heights.units <- rep("lines", n.row)
        heights.data <- as.list(1:n.row)

        heights.x[1] <- if (key$title == "") 0 else key$cex.title
        
        widths.x <- rep(key$between.column, n.col)
        widths.units <- rep("strwidth", n.col)
        widths.data <- as.list(rep("o", n.col))
        
        for (i in 1:column.blocks) {
            widths.x[(1:number.of.components-1)*3+1 +
                     (i-1)*3*number.of.components + i-1] <-
                         key$between/2
            
            widths.x[(1:number.of.components-1)*3+1 +
                     (i-1)*3*number.of.components + i+1] <-
                         key$between/2
        }
    
        
        for (i in 1:number.of.components) {

            cur <- components[[i]]

            id <- (1:column.blocks - 1) *
                (number.of.components*3 + 1) + i*3 - 1

            if (cur$type == "text") {
                which.name <- "W"
                for (ss in cur$pars$labels)
                    if (nchar(ss) > nchar(which.name)) which.name <- ss
                
                widths.x[id] <- max(cur$pars$cex)
                widths.units[id] <- "strwidth"

                for (idi in id) widths.data[[idi]] <- which.name

            }
            else if (cur$type == "rectangles") {
                widths.x[id] <- max(cur$pars$size)
            }
            else if (cur$type == "lines") {
                widths.x[id] <- max(cur$pars$size)
            }
            else if (cur$type == "points") {
                widths.x[id] <- max(cur$pars$cex)
            }
        }
        ## OK, layout set up, now to draw the key

        key.layout <- grid.layout(nrow = n.row, ncol = n.col,
                                  widths = unit(widths.x, widths.units, data=widths.data),
                                  heights = unit(heights.x, heights.units, data=heights.data),
                                  respect = TRUE)

        key.gf <- grid.frame(layout = key.layout, vp = vp,
                             gp = gpar(fontsize=9), # the value 9 used later
                             draw = FALSE)

        if (!key$transparent) {
            grid.place(key.gf,
                       grid.rect(gp=gpar(fill = key$background, col = key$border),
                                 draw = FALSE),
                       draw = FALSE, row = NULL, col = NULL)
        }
        else
            grid.place(key.gf,
                       grid.rect(gp=gpar(col=key$border), draw = FALSE),
                       draw = FALSE, row = NULL, col = NULL)

        ## Title
        if (!(key$title == ""))
            grid.place(key.gf, 
                       grid.text(label = key$title, draw = FALSE,
                                 gp = gpar(fontsize = 9 * key$cex.title)),  ## should do better than just 9
                       row=1, col = NULL, draw = FALSE)
        

        
        for (i in 1:number.of.components) {

            cur <- components[[i]]

            for (j in 1:cur$length) {

                colblck <- ceiling(j / rows.per.block)

                xx <- (colblck - 1) *
                    (number.of.components*3 + 1) + i*3 - 1

                yy <- j %% rows.per.block + 1
                if (yy == 1) yy <- rows.per.block + 1

                if (cur$type == "text") {
                    grid.place(key.gf, 
                               grid.text(label = cur$pars$labels[j],
                                         gp = gpar(col = cur$pars$col[j],
                                         font = cur$pars$font[j],
                                         fontsize = 9 * cur$pars$cex[j]),  ## should do better than just 9
                                         draw = FALSE),
                               row = yy, col = xx, draw = FALSE)
                    
                }
                else if (cur$type == "rectangles") {
                    grid.place(key.gf, 
                              grid.rect(width = cur$pars$size[j]/max(cur$pars$size),
                                        ## centred, unlike Trellis, due to aesthetic reasons !
                                        gp = gpar(fill = cur$pars$col[j], col = NULL), 
                                        draw = FALSE),
                              row = yy, col = xx, draw = FALSE)
                    
                    ## Need to make changes to support angle/density
                }
                else if (cur$type == "lines") {
                    if (cur$pars$type[j] == "l") {
                        grid.place(key.gf,
                                  grid.lines(x = c(0,1) * cur$pars$size[j]/max(cur$pars$size),
                                             ## ^^ this should be centered as well, but since the
                                             ## chances that someone would actually use this feature
                                             ## are astronomical, I'm leaving that for later.
                                             y = c(.5, .5),
                                             gp = gpar(col = cur$pars$col[j],
                                             lty = cur$pars$lty[j],
                                             lwd = cur$pars$lwd[j]),
                                             draw = FALSE),
                                  row = yy, col = xx, draw = FALSE)
                    }
                    else if (cur$pars$type[j] == "p") {
                        if (is.character(cur$pars$pch[j]))
                            grid.place(key.gf,
                                      grid.text(lab = cur$pars$pch[j], x = .5, y = .5,
                                                gp = gpar(col = cur$pars$col[j],
                                                fontsize = cur$pars$cex[j] * 10),
                                                draw = FALSE),
                                      row = yy, col = xx, draw = FALSE)
                        else {
                            grid.place(key.gf,
                                      grid.points(x=.5, y=.5, 
                                                  gp = gpar(col = cur$pars$col[j]),
                                                  size = unit(cur$pars$cex[j] * 2.5, "mm"),
                                                  pch = cur$pars$pch[j],
                                                  draw = FALSE),
                                      row = yy, col = xx, draw = FALSE)
                        }
                    }
                    else { # if (cur$pars$type[j] == "b" or "o") -- not differentiating
                        grid.place(key.gf, 
                                  grid.lines(x = c(0,1) * cur$pars$size[j]/max(cur$pars$size),
                                             ## ^^ this should be centered as well, but since the
                                             ## chances that someone would actually use this feature
                                             ## are astronomical, I'm leaving that for later.
                                             y = c(.5, .5),
                                             gp = gpar(col = cur$pars$col[j],
                                             lty = cur$pars$lty[j],
                                             lwd = cur$pars$lwd[j]),
                                             draw = FALSE),
                                  row = yy, col = xx, draw = FALSE)

                        if (is.character(cur$pars$pch[j]))
                            grid.place(key.gf, 
                                      grid.text(lab = cur$pars$pch[j],
                                                x = (1:key$divide-1)/(key$divide-1),
                                                y = rep(.5, key$divide),
                                                gp = gpar(col = cur$pars$col[j],
                                                fontsize = cur$pars$cex[j] * 10),
                                                draw = FALSE),
                                      row = yy, col = xx, draw = FALSE)
                        else
                            grid.place(key.gf, 
                                      grid.points(x = (1:key$divide-1)/(key$divide-1),
                                                  y = rep(.5, key$divide),
                                                  gp = gpar(col = cur$pars$col[j]),
                                                  size = unit(cur$pars$cex[j] * 2.5, "mm"),
                                                  pch = cur$pars$pch[j],
                                                  draw = FALSE),
                                      row = yy, col = xx, draw = FALSE)
                    }
                }
                else if (cur$type == "points") {
                    if (is.character(cur$pars$pch[j]))
                        grid.place(key.gf, 
                                  grid.text(lab = cur$pars$pch[j], x=.5, y=.5, 
                                            gp = gpar(col = cur$pars$col[j],
                                            fontsize = cur$pars$cex[j] * 10),
                                            draw = FALSE),
                                  row = yy, col = xx, draw = FALSE)
                    else {
                        grid.place(key.gf,
                                  grid.points(x=.5, y=.5, 
                                              gp = gpar(col = cur$pars$col[j]),
                                              size = unit(cur$pars$cex[j] * 2.5, "mm"),
                                              pch = cur$pars$pch[j],
                                              draw = FALSE),
                                  row = yy, col = xx, draw = FALSE)
                    }
                }

            }

        }

    }
    else stop("sorry, align=F not supported (yet ?)")


    if (draw)
        grid.draw(key.gf)

    key.gf
}


























draw.colorkey <- function(key, draw = FALSE, vp = NULL)
{
    if (!is.list(key)) stop("key must be a list")
    
    process.key <-
        function(col,
                 at,
                 tick.number = 7,
                 width = 2,
                 height = 1,
                 space = "right",
                 ...)
        {
            list(col = col,
                 at = at,
                 tick.number = tick.number,
                 width = width,
                 height = height,
                 space = space,
                 ...)
        }

    key <- do.call("process.key", key)
    
    ## Getting the locations/dimensions/centers of the rectangles
    key$at <- sort(key$at) ## should check if ordered
    if (length(key$at)!=length(key$col)+1) stop("length(col) must be length(at)-1")
    atrange <- range(key$at)
    scat <- .5 - key$height/2 + key$height *
        (key$at - atrange[1]) / diff(atrange)
    recnum <- length(scat)-1
    reccentre <- (scat[-1] + scat[-length(scat)]) / 2
    recdim <- diff(scat)

    cex <- 1
    col <- "black"
    font <- 1
    if (is.null(key$lab)) {
        at <- lpretty(atrange, key$tick.number)
        at <- at[at>=atrange[1] & at<=atrange[2]]
        labels <- as.character(at)
    }
    else if (is.character(key$lab) && length(key$lab)==length(key$at)) {
        at <- key$at
        labels <- as.character(key$lab)
    }
    else if (is.list(key$lab)) {
        if (!is.null(key$lab$at)) at <- key$lab$at
        if (!is.null(key$lab$lab)) labels <- as.character(key$lab$lab)
        if (!is.null(key$lab$cex)) cex <- key$lab$cex
        if (!is.null(key$lab$col)) col <- key$lab$col
        if (!is.null(key$lab$font)) font <- key$lab$font
    }
    else stop("malformed colorkey")
    labscat <- .5 - key$height/2 + key$height *
        (at - atrange[1]) / diff(atrange) # scales tick positions

    which.name <- "W"
    for (ss in labels)
        if (nchar(ss) > nchar(which.name)) which.name <- ss


    ## the tick marks for left and right should be modified

    
    if (key$space == "right") {

        widths.x <- c(key$width,1)
        widths.units <- rep("strwidth", 2)
        widths.data <- as.list(c("W", paste("--",which.name)))
        
        key.layout <-
            grid.layout(nrow = 1, ncol = 2,
                        widths = unit(widths.x, widths.units, data=widths.data))
        
        key.gf <- grid.frame(layout = key.layout, vp = vp,
                             gp = gpar(fontsize=8), #why 8 ?
                             draw = FALSE)
        
        for (c in seq(along=key$col))
            grid.pack(frame = key.gf, row = 1, col = 1,
                      grob =
                      grid.rect(y = unit(reccentre[c], "native"),
                                height = unit(recdim[c], "native"),
                                gp=gpar(fill=key$col[c],  col = NULL), draw = FALSE),
                      draw = FALSE)

        grid.pack(frame = key.gf, col = 1,
                  grob =
                  grid.rect(height = key$height,
                            gp=gpar(col="black"), draw = FALSE),
                  draw = FALSE)
        
        grid.pack(frame = key.gf, col = 1,
                  grob =
                  grid.text(label = paste("-", labels, sep = ""),
                            x = rep(1, length(labscat)),
                            y = labscat,
                            just = c("left","center"),
                            gp=gpar(col = col, fontsize = cex * 8, font = font),
                            draw = FALSE),
                  draw = FALSE)
    }
    else if (key$space == "left") {

        widths.x <- c(1,key$width)
        widths.units <- rep("strwidth", 2)
        widths.data <- as.list(c(paste("--",which.name), "W"))
        
        key.layout <-
            grid.layout(nrow = 1, ncol = 2,
                        widths = unit(widths.x, widths.units, data=widths.data))
        
        key.gf <- grid.frame(layout = key.layout, vp = vp,
                             gp = gpar(fontsize=8), #why 8 ?
                             draw = FALSE)
        
        for (c in seq(along=key$col))
            grid.pack(frame = key.gf, row = 1, col = 2,
                      grob =
                      grid.rect(y = unit(reccentre[c], "native"),
                                height = unit(recdim[c], "native"),
                                gp=gpar(fill=key$col[c],  col = NULL), draw = FALSE),
                      draw = FALSE)

        grid.pack(frame = key.gf, col = 2,
                  grob =
                  grid.rect(height = key$height,
                            gp=gpar(col="black"), draw = FALSE),
                  draw = FALSE)
        
        grid.pack(frame = key.gf, col = 2,
                  grob =
                  grid.text(label = paste(labels, "-", sep = ""),
                            x = rep(0, length(labscat)),
                            y = labscat,
                            just = c("right","center"),
                            gp=gpar(col = col, fontsize = cex * 8, font = font),
                            draw = FALSE),
                  draw = FALSE)
    }
    else if (key$space == "top") {

        heights.x <- c(1, .2, key$width)
        heights.units <- c("lines", "lines", "strwidth")
        heights.data <- as.list(c("a", "a", "W"))
        
        key.layout <-
            grid.layout(nrow = 3, ncol = 1,
                        heights = unit(heights.x, heights.units, data=heights.data))
        
        key.gf <- grid.frame(layout = key.layout, vp = vp,
                             gp = gpar(fontsize=8), #why 8 ?
                             draw = FALSE)
        
        for (c in seq(along=key$col))
            grid.pack(frame = key.gf, row = 3, col = 1,
                      grob =
                      grid.rect(x = unit(reccentre[c], "native"),
                                width = unit(recdim[c], "native"),
                                gp=gpar(fill=key$col[c],  col = NULL), draw = FALSE),
                      draw = FALSE)

        grid.pack(frame = key.gf, row = 3,
                  grob =
                  grid.rect(width = key$height,
                            gp=gpar(col="black"), draw = FALSE),
                  draw = FALSE)

        for (c in seq(along = labscat))
            grid.pack(frame = key.gf, row = 2,
                      grob =
                      grid.lines(x = unit(rep(labscat[c], 2), "native"),
                                 gp=gpar(col = col),
                                 draw = FALSE),
                                 draw = FALSE)

        grid.pack(frame = key.gf, row = 1,
                  grob =
                  grid.text(label = labels,
                            x = labscat,
                            just = c("centre","center"),
                            gp=gpar(col = col, fontsize = cex * 8, font = font),
                            draw = FALSE),
                  draw = FALSE)
    }
    else if (key$space == "bottom") {

        heights.x <- c(key$width, .2, 1)
        heights.units <- c("strwidth", "lines", "lines")
        heights.data <- as.list(c("W", "a", "a"))
        
        key.layout <-
            grid.layout(nrow = 3, ncol = 1,
                        heights = unit(heights.x, heights.units, data=heights.data))
        
        key.gf <- grid.frame(layout = key.layout, vp = vp,
                             gp = gpar(fontsize=8), #why 8 ?
                             draw = FALSE)
        
        for (c in seq(along=key$col))
            grid.pack(frame = key.gf, row = 1, col = 1,
                      grob =
                      grid.rect(x = unit(reccentre[c], "native"),
                                width = unit(recdim[c], "native"),
                                gp = gpar(fill=key$col[c], col = NULL), draw = FALSE),
                      draw = FALSE)

        grid.pack(frame = key.gf, row = 1,
                  grob =
                  grid.rect(width = key$height,
                            gp = gpar(col="black"), draw = FALSE),
                  draw = FALSE)

        for (c in seq(along = labscat))
            grid.pack(frame = key.gf, row = 2,
                      grob =
                      grid.lines(x = unit(rep(labscat[c], 2), "native"),
                                 gp=gpar(col = col),
                                 draw = FALSE),
                                 draw = FALSE)

        grid.pack(frame = key.gf, row = 3,
                  grob =
                  grid.text(label = labels,
                            x = labscat,
                            just = c("centre","center"),
                            gp=gpar(col = col, fontsize = cex * 8, font = font),
                            draw = FALSE),
                  draw = FALSE)
    }
    





    if (draw)
        grid.draw(key.gf)
    
    key.gf
}




























print.trellis <-
    function(x, position, split, more = FALSE, ...)
{
    if (is.null(dev.list())) trellis.device()
    else if (is.null(trellis.par.get())) trellis.device(device = .Device, new = FALSE)
    bg = trellis.par.get("background")$col
    new <- TRUE
    if (.lattice.print.more) new <- FALSE
    .lattice.print.more <<- more
    usual  <- (missing(position) & missing(split))
    ##if (!new && usual)
    ##    warning("more is relevant only when split/position is specified")
    
    if (!missing(position)) {
        if (length(position)!=4) stop("Incorrect value of position")
        if (new) {
            grid.newpage()
            grid.rect(gp = gpar(fill = bg, col = "transparent"))
        }
        push.viewport(viewport(x=position[1], y=position[2],
                               width=position[3]-position[1],
                               height=position[4]-position[2],
                               just=c("left","bottom")))
        
        if (!missing(split)) {
            if (length(split)!=4) stop("Incorrect value of split")

            push.viewport(viewport(layout = grid.layout(nrow=split[4],
                                   ncol = split[3])))
            push.viewport(viewport(layout.pos.row = split[2],
                                   layout.pos.col = split[1]))
        }
    }
    
    
    else if (!missing(split)) {
        
        if (length(split)!=4) stop("Incorrect value of split")
        if (new) {
            grid.newpage()
            grid.rect(gp = gpar(fill = bg, col = "transparent"))
        }
        push.viewport(viewport(layout = grid.layout(nrow=split[4],
                               ncol = split[3])))
        push.viewport(viewport(layout.pos.row = split[2],
                               layout.pos.col = split[1]))
    }
    
    panel <- # shall use "panel" in do.call
        if (is.function(x$panel)) x$panel 
        else if (is.character(x$panel)) get(x$panel)
        else eval(x$panel)

    x$strip <- 
        if (is.function(x$strip)) x$strip 
        else if (is.character(x$strip)) get(x$strip)
        else eval(x$strip)

    axis.line <- trellis.par.get("axis.line")
    number.of.cond <- length(x$condlevels)
    
    panel.width <- 1
    panel.height <- x$aspect.ratio
    layout.respect <- !x$aspect.fill

    if (!is.null(x$key)) {
        key.gf <- draw.key(x$key)
        key.space <-
            if ("space" %in% names(x$key)) x$key$space
            else if ("x" %in% names(x$key) ||
                     "corner" %in% names(x$key)) "inside"
            else "top"
    }
    else if (!is.null(x$colorkey)) {
        key.gf <- draw.colorkey(x$colorkey)
        key.space <- 
            if ("space" %in% names(x$colorkey)) x$colorkey$space
            else "right"
    }

    xaxis.col <-
        if (is.logical(x$x.scales$col)) axis.line$col
        else x$x.scales$col
    xaxis.font <-
        if (is.logical(x$x.scales$font)) 1
        else x$x.scales$font
    xaxis.cex <-
        x$x.scales$cex * x$fontsize.small / x$fontsize.normal
    xaxis.rot <-
        if (is.logical(x$x.scales$rot)) 0
        else x$x.scales$rot
    yaxis.col <-
        if (is.logical(x$y.scales$col)) axis.line$col
        else x$y.scales$col
    yaxis.font <-
        if (is.logical(x$y.scales$font)) 1
        else x$y.scales$font
    yaxis.cex <-
        x$y.scales$cex * x$fontsize.small / x$fontsize.normal
    yaxis.rot <-
        if (!is.logical(x$y.scales$rot)) x$y.scales$rot
        else if (x$y.scales$relation != "same" && is.logical(x$y.scales$labels)) 90
        else 0

    strip.col.default.bg <-
        rep(trellis.par.get("strip.background")$col,length=number.of.cond)
    strip.col.default.fg <-
        rep(trellis.par.get("strip.shingle")$col,length=number.of.cond)


    cond.max.level <- integer(number.of.cond)
    for(i in 1:number.of.cond) {
        cond.max.level[i] <- length(x$condlevels[[i]])
    }

    if(x$layout[1]==0) { # using device dimensions to
        ddim <- par("din") # calculate default layout
        device.aspect <- ddim[2]/ddim[1]
        panel.aspect <- if(layout.respect) panel.height else 1

        plots.per.page <- x$layout[2]
        m <- round(sqrt(x$layout[2] * device.aspect/panel.aspect))
        n <- ceiling(x$layout[2]/m)
        x$layout[1] <- n
        x$layout[2] <- m
    }
    else plots.per.page <- x$layout[1] * x$layout[2] 



    cols.per.page <- x$layout[1]
    rows.per.page <- x$layout[2]
    number.of.pages <- x$layout[3]
        
    if(cols.per.page>1)
        x.between <- rep(x$x.between, length = cols.per.page-1)
    if(rows.per.page>1) 
        y.between <- rep(x$y.between, length = rows.per.page-1)
    
    x.alternating <- rep(x$x.scales$alternating, length = cols.per.page)
    y.alternating <- rep(x$y.scales$alternating, length = rows.per.page)
    x.relation.same <- x$x.scales$relation == "same"
    y.relation.same <- x$y.scales$relation == "same"

    xlog <- x$x.scales$log
    ylog <- x$y.scales$log
    if (is.logical(xlog) && xlog) xlog <- 10
    if (is.logical(ylog) && ylog) ylog <- 10
    have.xlog <- !is.logical(xlog) || xlog
    have.ylog <- !is.logical(ylog) || ylog
    xlogbase <-
        if (is.numeric(xlog)) xlog
        else exp(1)
    ylogbase <-
        if (is.numeric(ylog)) ylog
        else exp(1)
    xlogpaste <-
        if (have.xlog) paste(as.character(xlog), "^", sep = "")
        else ""
    ylogpaste <-
        if (have.ylog) paste(as.character(ylog), "^", sep = "")
        else ""


    if (!is.logical(x$x.scales$at)) {  # i.e., at explicitly specified 
        check.x.overlap <- FALSE
        if (is.logical(x$x.scales$labels))
            if (have.xlog) {
                x$x.scales$labels <- as.character(x$x.scales$at)
                x$x.scales$at <- log(x$x.scales$at, xlogbase)
            }
            else x$x.scales$labels <- as.character(x$x.scales$at)
        else x$x.scales$labels <- as.character(x$x.scales$labels)
    }
    else check.x.overlap <- TRUE

    if (!is.null(x$x.scales$abbr)) { # this is for abbreviating long labels
        if (x$x.scales$abbr) 
            x$x.scales$labels <- 
                if (is.null(x$x.scales$minlen)) abbreviate(x$x.scales$labels)
                else abbreviate(x$x.scales$labels, x$x.scales$minlen)
    }

    if (!is.logical(x$y.scales$at)) {  # i.e., at explicitly specified 
        check.y.overlap <- FALSE
        if (is.logical(x$y.scales$labels))
            if (have.ylog) {
                x$y.scales$labels <- as.character(x$y.scales$at)
                x$y.scales$at <- log(x$y.scales$at, ylogbase)
            }
            else x$y.scales$labels <- as.character(x$y.scales$at)
        else x$y.scales$labels <- as.character(x$y.scales$labels)
    }
    else check.y.overlap <- TRUE

    if (!is.null(x$y.scales$abbr)) { # this is for abbreviating long labels
        if (x$y.scales$abbr) 
            x$y.scales$labels <- 
                if (is.null(x$y.scales$minlen)) abbreviate(x$y.scales$labels)
                else abbreviate(x$y.scales$labels, x$y.scales$minlen)
    }


    if (x.relation.same && is.logical(x$x.scales$at)) {
        x$x.scales$at <-
            lpretty(x$x.limits,
                    n = x$x.scales$tick.number)
        if (is.logical(x$x.scales$labels))
            x$x.scales$labels <- paste(xlogpaste, as.character(x$x.scales$at), sep = "")
        x$x.scales$labels <- as.character(x$x.scales$labels)
    }
    if (y.relation.same && is.logical(x$y.scales$at)) {
        x$y.scales$at <-
            lpretty(x$y.limits,
                    n = x$y.scales$tick.number)
        if (is.logical(x$y.scales$labels))
            x$y.scales$labels <- paste(ylogpaste, as.character(x$y.scales$at), sep = "")
        x$y.scales$labels <- as.character(x$y.scales$labels)
    }

    have.main <- !(is.null(x$main$label) || (is.character(x$main$label) && x$main$label==""))
    have.sub <- !(is.null(x$sub$label)   || (is.character(x$sub$label) && x$sub$label==""))
    have.xlab <- !(is.null(x$xlab$label) || (is.character(x$xlab$label) && x$xlab$label==""))
    have.ylab <- !(is.null(x$ylab$label) || (is.character(x$ylab$label) && x$ylab$label==""))

    
    ## Shall calculate the per page layout now:
    

    n.row <- rows.per.page * (number.of.cond + 3) + (rows.per.page-1) + 11
    ##       ^^^^^^^^^^^      ^^^^^^^^^^^^^^^^       ^^^^^^^^^^^^^^^    ^^
    ##          panels         rows per panel           between     see below
    ##               (2 for axes/ticks when relation!=same)

    ## the 11 things are as follows (top to bottom)
    ## 1/2 line space at top
    ## main
    ## key
    ## tick labels
    ## ticks
    ##
    ##   actual panels
    ##
    ## ticks
    ## tick labels
    ## xlab
    ## key
    ## sub
    ## 1/2 line space at bottom

    n.col <- 3 * cols.per.page + (cols.per.page-1) + 9 # similar

    ## the 9 things are as follows (left to right)
    ## 1/2 line space at left
    ## key
    ## ylab
    ## tick labels
    ## ticks
    ##
    ##   actual panels
    ##
    ## ticks
    ## tick labels
    ## key
    ## 1/2 line space at right

    heights.x <- rep(1, n.row)
    heights.units <- rep("lines", n.row)
    heights.data <- as.list(1:n.row)
                                        # not used now, but maybe needed if
                                        # rot = 90 for y-labels is supported

    widths.x <- rep(1, n.col)
    widths.units <- rep("lines", n.col)
    widths.data <- as.list(1:n.col) # this is required because allocating widths 
                                        # for y-axis annotation is more complicated

    ## fine tuning heights:


    heights.x[number.of.cond + 6 + (1:rows.per.page - 1) * (number.of.cond+4)] <-
        panel.height # for the panels
    heights.units[number.of.cond + 6 + (1:rows.per.page - 1) * (number.of.cond+4)] <-
        "null" # for the panels

    heights.x[number.of.cond + 7 + (1:rows.per.page - 1) * (number.of.cond+4)] <- 0
    ## This is for the x-axis ticks just below each panel if relation!="same"
    heights.x[number.of.cond + 8 + (1:rows.per.page - 1) * (number.of.cond+4)] <- 0
    ## This is for the x-axis labels just below each panel if relation!="same"

    heights.x[4] <- 0
    heights.x[5] <- 0 # tick axes/labels
    heights.x[n.row-4] <- 0
    heights.x[n.row-5] <- 0



    if (rows.per.page > 1)
        heights.x[number.of.cond + 9 + (
                                        (if (x$as.table) 1:(rows.per.page-1)
                                        else (rows.per.page-1):1)
                                        - 1)*(number.of.cond+4)] <-
                                            y.between
    ## y-between


    heights.x[1] <- 0.5
    heights.x[2] <- if (have.main) 2 * x$main$cex else 0
    if (have.main) {
        heights.units[2] <-  "strheight"
        heights.data[[2]] <- x$main$lab
    }



    heights.x[n.row] <- 0.5
    heights.x[n.row-1] <- if (have.sub) 2 * x$sub$cex else 0
    if (have.sub) {
        heights.units[n.row-1] <-  "strheight"
        heights.data[[n.row-1]] <- x$sub$lab
    }

    heights.x[3] <- 0 # for the key
    heights.x[n.row-2] <- 0 # key




    if (x$x.scales$draw) {

        if (x.relation.same) {

            heights.x[5] <- x$x.scales$tck * 0.453
            heights.x[n.row-5] <- x$x.scales$tck * 0.453

            if (xaxis.rot == 0) {

                if (any(x.alternating==2)) 
                    heights.x[4] <- xaxis.cex

                if (any(x.alternating==1)) 
                    heights.x[n.row-4] <- xaxis.cex

            }
            else {

                which.name <- "W"
                for (ss in x$x.scales$labels)
                    if (nchar(ss) > nchar(which.name)) which.name <- ss
                if (any(x.alternating==2)) {
                    heights.x[4] <- xaxis.cex * abs(sin(xaxis.rot * pi /180))
                    heights.units[4] <- "strwidth"
                    heights.data[[4]] <- which.name
                }
                if (any(x.alternating==1)) {
                    heights.x[n.row-4] <- xaxis.cex * abs(sin(xaxis.rot * pi /180))
                    heights.units[n.row-4] <- "strwidth"
                    heights.data[[n.row-4]] <- which.name
                }
            }
        }
        else { # relation != same

            if (xaxis.rot == 0) {

                heights.x[number.of.cond + 7 + (1:rows.per.page - 1)*(number.of.cond+4)] <- x$x.scales$tck * 0.2
                heights.x[number.of.cond + 8 + (1:rows.per.page - 1)*(number.of.cond+4)] <- xaxis.cex

            }
            else {

                which.name <- "W"

                if (is.logical(x$x.scales$at)) {

                    for (sc in x$x.limits) {
                        at <- lpretty(sc, n = x$x.scales$tick.number)
                        labs <- paste(xlogpaste, as.character(at), sep = "")
                        for (ss in labs)
                            if (nchar(ss) > nchar(which.name)) which.name <- ss
                    }
                }
                else
                    for (ss in x$x.scales$labels)
                        if (nchar(ss) > nchar(which.name)) which.name <- ss
                

                heights.x[number.of.cond + 7 + (1:rows.per.page - 1)*(number.of.cond+4)] <- x$x.scales$tck * 0.3
                heights.x[number.of.cond + 8 + (1:rows.per.page - 1)*(number.of.cond+4)] <- xaxis.cex * abs(sin(xaxis.rot * pi /180))
                heights.units[number.of.cond + 8 + (1:rows.per.page - 1)*(number.of.cond+4)] <- "strwidth"
                heights.data[number.of.cond + 8 + (1:rows.per.page - 1)*(number.of.cond+4)] <- which.name

            }
        }
    }

    heights.x[n.row-3] <- if (have.xlab) 2 * x$xlab$cex else 0 # xlab
    if (have.xlab) {
        heights.units[n.row-3] <-  "strheight"
        heights.data[[n.row-3]] <- x$xlab$lab
    }

    ## this is if strip=F -- strips not to be drawn
    for(crr in 1:number.of.cond)
        heights.x[number.of.cond + 6 + (1:rows.per.page - 1)*(number.of.cond+4) - crr] <-
            if (is.logical(x$strip)) 0  # which means strip = F, strips not to be drawn
            else 1.1 * x$par.strip.text$cex * x$par.strip.text$lines

    ## fine tuning widths:
    ##----------------------------------------------------------------------------------

    widths.x[3] <- if (have.ylab) 2 * x$ylab$cex else 0 # ylab
    if (have.ylab) {
        widths.units[3] <-  "strheight"
        widths.data[[3]] <- x$ylab$lab
    }


    widths.x[(1:cols.per.page - 1)*4 + 8] <-
        panel.width # for the panels
    widths.units[(1:cols.per.page - 1)*4 + 8] <-
        "null" # for the panels


    widths.x[(1:cols.per.page - 1)*4 + 7] <- 0
    widths.x[(1:cols.per.page - 1)*4 + 6] <- 0
    ## For y-axes labels and ticks to the left of each panel when relation != "same"
    ## (might change later)

    widths.x[4] <- 0
    widths.x[5] <- 0 #ticks/labels
    widths.x[n.col-2] <- 0
    widths.x[n.col-3] <- 0

    if (cols.per.page > 1)
        widths.x[(1:(cols.per.page-1) - 1)*4 + 9] <- x.between
    ## x-between

    widths.x[1] <- 0.5
    widths.x[n.col] <- 0.5
    widths.x[2] <- 0 # key - left
    widths.x[n.col-1] <- 0 # key - right

    ## next part of the code decides how much space to leave for y-labels

    if (x$y.scales$draw) {

        if (y.relation.same) {

            widths.x[5] <- x$y.scales$tck * 0.453 ## tck = 1 will be 0.453 "lines"
            widths.x[n.col-3] <- x$y.scales$tck * 0.453

            if (abs(yaxis.rot) == 90) {

                if(any(y.alternating==1)) 
                    widths.x[4] <- yaxis.cex

                if (any(y.alternating==2)) 
                    widths.x[n.col-2] <- yaxis.cex

            }
            else {
                
                which.name <- "W"
                for (ss in x$y.scales$labels)
                    if (nchar(ss) > nchar(which.name)) which.name <- ss
                if(any(y.alternating==1)) {
                    widths.x[4] <- yaxis.cex * abs(cos(yaxis.rot * pi /180))
                    widths.units[4] <- "strwidth"
                    widths.data[[4]] <- which.name
                }
                if (any(y.alternating==2)) {
                    widths.x[n.col-2] <- yaxis.cex * abs(cos(yaxis.rot * pi /180))
                    widths.units[n.col-2] <- "strwidth"
                    widths.data[[n.col-2]] <- which.name
                }
            }
        }
        else { # relation != same
            if (abs(yaxis.rot) == 90) {

                widths.x[(1:cols.per.page - 1)*4 + 6] <- yaxis.cex
                widths.x[(1:cols.per.page - 1)*4 + 7] <- x$y.scales$tck * 0.2

            }
            else {

                which.name <- "W"

                if (is.logical(x$y.scales$at)) {

                    for (sc in x$y.limits) {
                        at <- lpretty(sc, n = x$y.scales$tick.number)
                        labs <- paste(ylogpaste, as.character(at), sep = "")
                        for (ss in labs)
                            if (nchar(ss) > nchar(which.name)) which.name <- ss
                    }
                }
                else
                    for (ss in x$y.scales$labels)
                        if (nchar(ss) > nchar(which.name)) which.name <- ss
                

                widths.x[(1:cols.per.page - 1)*4 + 7] <- x$y.scales$tck * 0.3
                widths.x[(1:cols.per.page - 1)*4 + 6] <- yaxis.cex * abs(cos(yaxis.rot * pi /180))
                widths.units[(1:cols.per.page - 1)*4 + 6] <- "strwidth"
                widths.data[(1:cols.per.page - 1)*4 + 6] <- which.name
                
            }
            
        }
        
    }


    if (!is.null(x$key) || !is.null(x$colorkey)) {
            
        if (key.space == "left") {
            widths.x[2] <- 1
            widths.units[2] <- "grobwidth"
            widths.data[[2]] <- key.gf
        }
        else if (key.space == "right") {
            widths.x[n.col-1] <- 1
            widths.units[n.col-1] <- "grobwidth"
            widths.data[[n.col-1]] <- key.gf
        }
        else if (key.space == "top") {
            heights.x[3] <- 1
            heights.units[3] <- "grobheight"
            heights.data[[3]] <- key.gf
        }
        else if (key.space == "bottom") {
            heights.x[n.row-2] <- 1
            heights.units[n.row-2] <- "grobheight"
            heights.data[[n.row-2]] <- key.gf
        }
        
    }
    

    ## Constructing the layout:


    page.layout <- grid.layout(nrow = n.row, ncol = n.col,
                               widths = unit(widths.x, widths.units, data=widths.data),
                               heights = unit(heights.x, heights.units, data=heights.data),
                               respect = layout.respect)

    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1
        
    for(page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0)) {
                
            if(usual) {
                if (new) grid.newpage()
                grid.rect(gp = gpar(fill = bg, col = "transparent"))
                new <- TRUE
            }

            push.viewport(viewport(layout = page.layout,
                                   gp = gpar(fontsize = x$fontsize.normal,
                                   col = axis.line$col,
                                   lty = axis.line$lty,
                                   lwd = axis.line$lwd)))

            if (have.main)
                grid.text(label = x$main$label,
                          gp = gpar(col = x$main$col, font = x$main$font, 
                          fontsize = x$fontsize.normal * x$main$cex),
                          vp = viewport(layout.pos.row = 2))
                    
                    
            if (have.sub)
                grid.text(label = x$sub$label,
                          gp = gpar(col = x$sub$col, font = x$sub$font, 
                          fontsize = x$fontsize.normal * x$sub$cex),
                          vp = viewport(layout.pos.row = n.row-1))
                    
                    
            if (have.xlab) 
                grid.text(label = x$xlab$label,
                          gp = gpar(col = x$xlab$col, font = x$xlab$font, 
                          fontsize = x$fontsize.normal * x$xlab$cex), 
                          vp = viewport(layout.pos.row = n.row - 3, layout.pos.col = c(6, n.col - 4)))
                
            if (have.ylab)
                grid.text(label = x$ylab$label, rot = 90,
                          gp = gpar(col = x$ylab$col, font = x$ylab$font, 
                          fontsize = x$fontsize.normal * x$ylab$cex),
                          vp = viewport(layout.pos.col = 3, layout.pos.row = c(6, n.row - 6)))
            
            for (row in 1:rows.per.page)
                for (column in 1:cols.per.page)

                    if (!any(cond.max.level-cond.current.level<0) &&
                        (row-1) * cols.per.page + column <= plots.per.page) {

                        if (!is.list(x$panel.args[[panel.number]]))
                            ## corr to skip = T or extra plots
                            panel.number <- panel.number + 1
                            
                        else {
                                
                            actual.row <- if (x$as.table)
                                (rows.per.page-row+1) else row
                            ## this gives the row position from the bottom


                            pos.row <- 6 + number.of.cond + 
                                (rows.per.page - actual.row) *
                                (number.of.cond + 4)
                            pos.col <- (column-1) * 4 + 8

                            xscale <-
                                if(x.relation.same)
                                    x$x.limits
                                else x$x.limits[[panel.number]]
                            yscale <- 
                                if(y.relation.same)
                                    x$y.limits
                                else x$y.limits[[panel.number]]
                            
                                
                            push.viewport(viewport(layout.pos.row = pos.row,
                                                   layout.pos.col = pos.col,
                                                   xscale = xscale,
                                                   yscale = yscale,
                                                   clip = TRUE,
                                                   gp = gpar(fontsize =
                                                   x$fontsize.normal)))


                            pargs <- c(x$panel.args[[panel.number]],
                                       x$panel.args.common)
                            if (!("..." %in% names(formals(panel))))
                                pargs <- pargs[names(formals(panel))]
                            do.call("panel", pargs)

                            grid.rect()

                            pop.viewport()

                            ## next few lines deal with drawing axes
                            ## as appropriate
                            ## when relation != same, axes drawn for
                            ## each panel:
                            
                            ## X-axis
                            if (!x.relation.same && x$x.scales$draw) {

                                axs <- x$x.scales

                                if (is.logical(axs$at)) {
                                    axs$at <- lpretty(xscale, n = axs$tick.number)
                                    axs$labels <- paste(xlogpaste, as.character(axs$at), sep = "")
                                }

                                ok <- (axs$at>=xscale[1] & axs$at<=xscale[2])

                                push.viewport(viewport(layout.pos.row = pos.row+1,
                                                       layout.pos.col = pos.col,
                                                       xscale = xscale))


                                for(tt in seq(along=axs$at)[ok])
                                    grid.lines(x = unit(rep(axs$at[tt],2), "native"),
                                               y = unit(c(0.1,1), "npc"),
                                               gp = gpar(col = xaxis.col))
                                pop.viewport()


                                grid.text(label = axs$label[ok],
                                          x = unit(axs$at[ok], "native"),
                                          just = "centre",
                                          ## just = if (xaxis.rot == 0) c("centre", "top")
                                          ## else c("right", "centre"),
                                          rot = xaxis.rot,
                                          check.overlap = check.x.overlap,
                                          gp = gpar(col = xaxis.col, font = xaxis.font, 
                                          fontsize = axs$cex * x$fontsize.small),
                                          vp = viewport(layout.pos.row = pos.row+2,
                                          layout.pos.col = pos.col, xscale = xscale))


                            }
                            ## Y-axis
                            if (!y.relation.same && x$y.scales$draw) {

                                axs <- x$y.scales

                                if (is.logical(axs$at)) {
                                    axs$at <- lpretty(yscale, n = axs$tick.number)
                                    axs$labels <- paste(ylogpaste, as.character(axs$at), sep = "")
                                }

                                ok <- (axs$at>=yscale[1] & axs$at<=yscale[2])

                                push.viewport(viewport(layout.pos.row = pos.row,
                                                       layout.pos.col = pos.col-1,
                                                       yscale = yscale))


                                for(tt in seq(along=axs$at)[ok])
                                    grid.lines(y = unit(rep(axs$at[tt],2), "native"),
                                               x = unit(c(0.1,1), "npc"),
                                               gp = gpar(col = yaxis.col))
                                pop.viewport()


                                grid.text(label = axs$label[ok],
                                          y = unit(axs$at[ok], "native"),
                                          just = "centre",
                                          ## just = if (xaxis.rot == 0) c("centre", "top")
                                          ## else c("right", "centre"),
                                          rot = yaxis.rot,
                                          check.overlap = check.y.overlap,
                                          gp = gpar(col = yaxis.col, font = yaxis.font, 
                                          fontsize = axs$cex * x$fontsize.small),
                                          vp = viewport(layout.pos.row = pos.row,
                                          layout.pos.col = pos.col-2, yscale = yscale))


                            }

                            ## When relation = same, axes drawn based on value of alternating
                            if (y.relation.same && x$y.scales$draw) {
                                
                                ## Y-axis to the left
                                if (column == 1) {
                                    axs <- x$y.scales

                                    ok <- (axs$at>=yscale[1] & axs$at<=yscale[2])

                                    push.viewport(viewport(layout.pos.row = pos.row,
                                                           layout.pos.col = pos.col-3,
                                                           yscale = yscale))


                                    for(tt in seq(along=axs$at)[ok])
                                        grid.lines(y = unit(rep(axs$at[tt],2), "native"),
                                                   x = unit(c(0.5,1), "npc"),
                                                   gp = gpar(col = yaxis.col))
                                    pop.viewport()

                                    if (y.alternating[actual.row]==1) {
                                        ##grid.rect(gp = gpar(fill = "cyan"),
                                        ##          vp = viewport(layout.pos.row = pos.row,
                                        ##          layout.pos.col = pos.col-4, yscale = yscale))
                                        grid.text(label = axs$label[ok],
                                                  x = if (abs(yaxis.rot) != 90) 1 else .5,
                                                  y = unit(axs$at[ok], "native"),
                                                  just = if (yaxis.rot != 90) c("right", "centre")
                                                  else c("centre", "centre"),
                                                  rot = yaxis.rot,
                                                  check.overlap = check.y.overlap,
                                                  gp = gpar(col = yaxis.col, font = yaxis.font, 
                                                  fontsize = axs$cex * x$fontsize.small),
                                                  vp = viewport(layout.pos.row = pos.row,
                                                  layout.pos.col = pos.col-4, yscale = yscale))
                                    }
                                }
                                
                                ## Y-axis to the right
                                if (column == cols.per.page) {

                                    axs <- x$y.scales

                                    ok <- (axs$at>=yscale[1] & axs$at<=yscale[2])

                                    push.viewport(viewport(layout.pos.row = pos.row,
                                                           layout.pos.col = pos.col+1,
                                                           yscale = yscale))


                                    for(tt in seq(along=axs$at)[ok])
                                        grid.lines(y = unit(rep(axs$at[tt],2), "native"),
                                                   x = unit(c(0,0.5), "npc"),
                                                   gp = gpar(col = yaxis.col))
                                    pop.viewport()

                                    if (y.alternating[actual.row]==2) {
                                        ##grid.rect(gp = gpar(fill = "cyan"),
                                        ##          vp = viewport(layout.pos.row = pos.row,
                                        ##          layout.pos.col = pos.col+2, yscale = yscale))
                                        grid.text(label = axs$label[ok],
                                                  x = if (abs(yaxis.rot) != 90) 0 else .5,
                                                  y = unit(axs$at[ok], "native"),
                                                  just = if (yaxis.rot != 90) c("left", "centre")
                                                  else c("centre", "centre"),
                                                  rot = yaxis.rot,
                                                  check.overlap = check.y.overlap,
                                                  gp = gpar(col = yaxis.col, font = yaxis.font, 
                                                  fontsize = axs$cex * x$fontsize.small),
                                                  vp = viewport(layout.pos.row = pos.row,
                                                  layout.pos.col = pos.col+2, yscale = yscale))
                                    }
                                }
                                
                            }
                                    
                                    
                                
                            ## X-axis to the bottom
                            if (x.relation.same && x$x.scales$draw) {

                                if (actual.row == 1) {

                                    axs <- x$x.scales

                                    ok <- (axs$at>=xscale[1] & axs$at<=xscale[2])

                                    push.viewport(viewport(layout.pos.row = pos.row+3,
                                                           layout.pos.col = pos.col,
                                                           xscale = xscale))



                                    for(tt in seq(along=axs$at)[ok])
                                        grid.lines(x = unit(rep(axs$at[tt],2), "native"),
                                                   y = unit(c(0.5,1), "npc"),
                                                   gp = gpar(col = xaxis.col))
                                    pop.viewport()

                                    if (x.alternating[column]==1) {
                                        ##grid.rect(gp = gpar(fill = "cyan"),
                                        ##          vp = viewport(layout.pos.row = pos.row + 4,
                                        ##          layout.pos.col = pos.col, xscale = xscale))
                                        grid.text(label = axs$label[ok],
                                                  y = if (xaxis.rot != 0) 1 else .5,
                                                  x = unit(axs$at[ok], "native"),
                                                  just = if (xaxis.rot != 0) c("right", "centre")
                                                  else c("centre", "centre"),
                                                  rot = xaxis.rot,
                                                  check.overlap = check.x.overlap,
                                                  gp = gpar(col = xaxis.col, font = xaxis.font, 
                                                  fontsize = axs$cex * x$fontsize.small),
                                                  vp = viewport(layout.pos.row = pos.row + 4,
                                                  layout.pos.col = pos.col, xscale = xscale))
                                    }
                                }
                            }
                                    
                            ##-------------------------


                            if (!is.logical(x$strip)) # logical ==> FALSE
                                for(i in 1:number.of.cond)
                                {
                                    push.viewport(viewport(layout.pos.row = pos.row-i,
                                                           layout.pos.col = pos.col,
                                                           clip = TRUE,
                                                           gp = gpar(fontsize = x$fontsize.normal)))
                                    
                                    grid.rect()
                                    x$strip(which.given = i,
                                            which.panel = cond.current.level,
                                            var.name = names(x$condlevels),
                                            factor.levels = if (!is.list(x$cond[[i]]))
                                            x$cond[[i]] else NULL,
                                            shingle.intervals = if (is.list(x$cond[[i]]))
                                            do.call("rbind", x$cond[[i]]) else NULL,
                                            ##x = x$condlevel[[i]],
                                            ##level = cond.current.level[i],
                                            ##name = names(x$cond)[i],
                                            bg = strip.col.default.bg[i],
                                            fg = strip.col.default.fg[i],
                                            par.strip.text = x$par.strip.text)
                                    
                                    pop.viewport()
                                            
                                }
                            
                            
                            ## X-axis at top
                            if (x.relation.same && x$x.scales$draw)
                                if (actual.row == rows.per.page) {

                                    axs <- x$x.scales

                                    ok <- (axs$at>=xscale[1] & axs$at<=xscale[2])

                                    push.viewport(viewport(layout.pos.row = pos.row - 1 - 
                                                           number.of.cond,
                                                           layout.pos.col = pos.col,
                                                           xscale = xscale))


                                    for(tt in seq(along=axs$at)[ok])
                                        grid.lines(x = unit(rep(axs$at[tt],2), "native"),
                                                   y = unit(c(0,0.5), "npc"),
                                                   gp = gpar(col = xaxis.col))
                                    pop.viewport()

                                    if (x.alternating[column]==2)
                                        grid.text(label = axs$label[ok],
                                                  y = if (xaxis.rot != 0) 0 else .5,
                                                  x = unit(axs$at[ok], "native"),
                                                  just = if (xaxis.rot != 0) c("left", "centre")
                                                  else c("centre", "centre"),
                                                  rot = xaxis.rot,
                                                  check.overlap = check.x.overlap,
                                                  gp = gpar(col = xaxis.col, font = xaxis.font, 
                                                  fontsize = axs$cex * x$fontsize.small),
                                                  vp = viewport(layout.pos.row = pos.row - 2 - 
                                                  number.of.cond,
                                                  layout.pos.col = pos.col, xscale = xscale))

                                }

                                
                            cond.current.level <- cupdate(cond.current.level,
                                                          cond.max.level)
                            panel.number <- panel.number + 1
                            
                        }
                        
                    }
            
            
            if (!is.null(x$key) || !is.null(x$colorkey)) {
                
                if (key.space == "left") {
                    push.viewport(viewport(layout.pos.col = 2,
                                  layout.pos.row = c(6, n.row-6)))
                    grid.draw(key.gf)
                    pop.viewport()
                }
                else if (key.space == "right") {
                    push.viewport(viewport(layout.pos.col = n.col-1,
                                  layout.pos.row = c(6, n.row-6)))
                    grid.draw(key.gf)
                    pop.viewport()
                    }
                else if (key.space == "top") {
                    push.viewport(viewport(layout.pos.row = 3,
                                           layout.pos.col = c(6,n.col-4)))
                    grid.draw(key.gf)
                    pop.viewport()
                }
                else if (key.space == "bottom") {
                    push.viewport(viewport(layout.pos.row = n.row - 2,
                                           layout.pos.col = c(6,n.col-4)))
                    grid.draw(key.gf)
                    pop.viewport()
                }
                else if (key.space == "inside") {
                    
                    push.viewport(viewport(layout.pos.row = c(1, n.row),
                                           layout.pos.col = c(1, n.col)))
                    
                    if (is.null(x$key$corner)) x$key$corner <- c(0,1)
                    if (is.null(x$key$x)) x$key$x <- x$key$corner[1]
                    if (is.null(x$key$y)) x$key$y <- x$key$corner[2]
                    
                    if (all(x$key$corner == c(0,1))) {
                        
                        push.viewport(viewport(layout = grid.layout(nrow = 3, ncol = 3,
                                               widths = unit(c(x$key$x, 1, 1),
                                               c("npc", "grobwidth", "null"),
                                               list(1, key.gf, 1)),
                                               heights = unit(c(1-x$key$y, 1, 1),
                                               c("npc", "grobheight", "null"),
                                               list(1, key.gf, 1)))))
                        
                        push.viewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
                        
                        grid.draw(key.gf)
                        
                        pop.viewport()
                        pop.viewport()
                        
                    }
                    
                    
                    if (all(x$key$corner == c(1,1))) {
                        
                        push.viewport(viewport(layout = grid.layout(nrow = 3, ncol = 3,
                                               heights = unit(c(1-x$key$y, 1, 1),
                                               c("npc", "grobheight", "null"),
                                               list(1, key.gf, 1)),
                                               widths = unit(c(1, 1, 1-x$key$x),
                                               c("null", "grobwidth", "npc"),
                                               list(1, key.gf, 1)))))
                        
                        push.viewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
                        
                        grid.draw(key.gf)

                        pop.viewport()
                        pop.viewport()
                        
                    }
                    

                    if (all(x$key$corner == c(0,0))) {
                    
                        push.viewport(viewport(layout = grid.layout(nrow = 3, ncol = 3,
                                               widths = unit(c(x$key$x, 1, 1),
                                               c("npc", "grobwidth", "null"),
                                               list(1, key.gf, 1)),
                                               heights = unit(c(1,1,x$key$y),
                                               c("null", "grobheight", "npc"),
                                               list(1, key.gf, 1)))))
                        
                        push.viewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
                        
                        grid.draw(key.gf)
                        
                        pop.viewport()
                        pop.viewport()
                        
                    }
                    
                    
                    if (all(x$key$corner == c(1,0))) {
                            
                        push.viewport(viewport(layout=grid.layout(nrow = 3, ncol = 3,
                                               widths = unit(c(1, 1, 1-x$key$x),
                                               c("null", "grobwidth", "npc"),
                                               list(1, key.gf, 1)),
                                               heights = unit(c(1, 1, x$key$y),
                                               c("null", "grobheight", "npc"),
                                               list(1, key.gf, 1)))))
                        
                        push.viewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
                        
                        grid.draw(key.gf)
                        
                        pop.viewport()
                        pop.viewport()
                        
                    }
                    
                    
                    
                    pop.viewport()
                    
                }
                
            }
            
            push.viewport(viewport(layout.pos.row = c(1, n.row),
                                   layout.pos.col = c(1, n.col)))
            if(!is.null(x$page)) x$page(page.number)                
            pop.viewport()
            
            pop.viewport()
            
            
        }

    if (!missing(position)) {
        if (!missing(split)) {
            pop.viewport()
            pop.viewport()
        }
        pop.viewport()
    }
    else if (!missing(split)) {
        pop.viewport()
        pop.viewport()
    }
    invisible()
}



### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



prepanel.default.qq <-
    function(x, y, ...)
{
    list(xlim = range(x, y),
         ylim = range(x, y),
         dx = 1,
         dy = 1)
}




panel.qq <-
    function(...)
{
    reference.line <- trellis.par.get("reference.line")
    panel.abline(0,1,
                 col = reference.line$col,
                 lty = reference.line$lty,
                 lwd = reference.line$lwd)
    panel.xyplot(...)

}



qq <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.qq",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             f.value = ppoints,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ## dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing
    
    form <- latticeParseFormula(formula, data)
    cond <- form$condition
    number.of.cond <- length(cond)
    y <- form$left
    x <- form$right
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    if(subscripts) subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    y <- y[subset, drop = TRUE]
    if (subscripts) subscr <- subscr[subset, drop = TRUE]
    
    x <- as.numeric(x)
    y <- as.factorOrShingle(y)
    is.f.y <- is.factor(y)
    num.l.y <- nlevels(y)
    if (num.l.y!=2) stop("y must have exactly 2 levels")

    if(missing(xlab)) xlab <-
        if (is.f.y) unique(levels(y))[1]
        else paste("y:", as.character(unique(levels(y)[[1]])))
    
    if(missing(ylab)) ylab <-
        if (is.f.y) unique(levels(y))[y]
        else paste("y:", as.character(unique(levels(y)[[2]])))


    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <-
            if (is.f.y) unique(levels(y))[1]
            else paste("y:", as.character(unique(levels(y)[[1]])))
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <-
            if (is.f.y) unique(levels(y))[y]
            else paste("y:", as.character(unique(levels(y)[[2]])))

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo, 
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used: completed later

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        ## x <- log(x, xbase)  later, in panel.args
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        ylog <- foo$y.scales$log
        ybase <-
            if (is.logical(ylog)) 10
            else if (is.numeric(ylog)) ylog
            else if (ylog == "e") exp(1)

        ## y <- log(y, ybase)
        if (have.ylim) ylim <- log(ylim, ybase)
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)|is.na(y)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args


    foo$panel.args.common <- dots
    if (subscripts) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    if (any(id)) {

                        if (is.f.y) {
                            tx <- x[id]
                            ty <- as.numeric(y[id])
                            x.val <- tx[ty==1]
                            y.val <- tx[ty==2]
                        }
                        else {
                            tx <- x[id]
                            ty <- y[id]
                            ly <- levels(y)
                            x.val <- tx[ty>=ly[[1]][1] & ty <=ly[[1]][2]]
                            y.val <- tx[ty>=ly[[2]][1] & ty <=ly[[2]][2]]
                        }
                        n <- max(length(x.val), length(y.val))
                        p  <- f.value(n)
                        foo$panel.args[[panel.number]] <-
                            list(x = quantile(x = x.val, probs = p),
                                 y = quantile(x = y.val, probs = p))

                    }
                    else
                        foo$panel.args[[panel.number]] <-
                            list(x = numeric(0), y = numeric(0))

                    if (subscripts)
                        foo$panel.args[[panel.number]]$subscripts <-
                            subscr[id]

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.qq,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}




### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA




panel.qqmathline <-
    function(y, distribution, ...)
{
    if (length(y) > 0) {
        yy <- quantile(y, c(0.25, 0.75))
        xx <- distribution(c(0.25, 0.75))
        r <- diff(yy)/diff(xx)
        panel.abline(c( yy[1]-xx[1]*r , r), ...)
    }
}


prepanel.qqmathline <-
    function(y, distribution, f.value = ppoints, ...)
{
    yy <- quantile(y, c(0.25, 0.75))
    xx <- distribution(c(0.25, 0.75))
    n <- length(y)
    list(ylim = range(y), xlim = range(distribution(f.value(n))),
         dx = diff(xx), dy = diff(yy))
}


prepanel.default.qqmath <-
    function(...)
    prepanel.default.xyplot(...)


panel.qqmath <-
    function(...) panel.xyplot(...)


qqmath <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.qqmath",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             f.value = ppoints,
             distribution = qnorm,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ## dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    distribution.name <- deparse(substitute(distribution))
    distribution <-
        if (is.function(distribution)) distribution 
        else if (is.character(distribution)) get(distribution)
        else eval(distribution)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing
    
    form <- latticeParseFormula(formula, data)
    cond <- form$condition
    number.of.cond <- length(cond)
    x <- form$right
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    if(subscripts) subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    if (subscripts) subscr <- subscr[subset, drop = TRUE]

    if(missing(ylab)) ylab <- form$right.name
    if(missing(xlab)) xlab <- distribution.name
    if (is.shingle(x)) stop("x cannot be a shingle")
    x <- as.numeric(x)

    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- deparse(substitute(distribution))
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <- form$right.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo, 
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        warning("cannot use log scale for y, use different distribution")
        foo$y.scales$log <- FALSE
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    foo$panel.args.common <- dots
    if (subscripts) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }
                    if (any(id)) {
                        foo$panel.args[[panel.number]] <-
                            list(x = distribution(f.value(length(x[id]))), 
                                 y = quantile(x[id], f.value(length(x[id]))))
                    }
                    else
                        foo$panel.args[[panel.number]] <-
                            list(x = numeric(0), y = numeric(0))

                    if (subscripts)
                        foo$panel.args[[panel.number]]$subscripts <-
                            subscr[id]

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.qqmath,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots,
                               distribution = distribution))

    class(foo) <- "trellis"
    foo
}







rfs <-
    function(model, layout = c(2,1), xlab = "f-value", ylab = NULL,
             distribution = qunif,
             panel = function(...) {panel.grid(); panel.qqmath(...)},
             prepanel = NULL, strip = TRUE, ...)
{

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    fitval <- fitted.values(model) - mean(fitted.values(model))
    resids <- residuals(model)
    
    nf <- length(fitval)
    nr <- length(resids)
    
    data <- list(y = c( fitval, resids),
                 f = c( rep("Fitted Values minus Mean", nf),
                 rep("Residuals", nr)))

    qqmath(~y|f, data = data, layout = layout, xlab = xlab, ylab = ylab,
           distribution = distribution, panel = panel,
           prepanel = prepanel, strip = strip, ...)
}


### Copyright 2001-2002 Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA


lattice.theme <- list()

col.whitebg <- function()
    list(background = list(col="transparent"),
         bar.fill = list(col="#c8ffc8"),
         box.rectangle = list(col="darkgreen"),
         box.umbrella = list(col="darkgreen"),
         dot.line = list(col="#e8e8e8"),
         dot.symbol = list(col="darkgreen"),
         plot.line = list(col="darkgreen"),
         plot.symbol = list(col="darkgreen"),
         ##regions=list(col=rev(hsv(h=250:349/1000, v=30:129/150,s=.5,
         ##gamma = .6)))
         regions = list(col = heat.colors(100)),
         strip.background = list(col = c("#ff7f00", "#00ff00", "#00ffff",
                                 "#0080ff", "#ff00ff", "#ff0000", "#ffff00")),
         strip.background = list(col = c("#ffe5cc", "#ccffcc", "#ccffff",
                                 "#cce6ff", "#ffccff", "#ffcccc", "#ffffcc")),
         reference.line = list(col="#e8e8e8"),
         superpose.line = list(col = c("darkgreen","red","royalblue",
                               "brown","orange","turquoise", "orchid"),
         lty = 1:7),
         superpose.symbol = list(pch = c(1,3,6,0,5,16,17), cex = rep(.7, 7),
         col = c("darkgreen","red","royalblue",
                                 "brown","orange","turquoise", "orchid")))


canonical.theme <- function(name = "null device", color = TRUE)
{
    ## For the purpose of this function, the only differences in the
    ## settings/themes arise from the difference in the default
    ## colors. So, I will first set up the appropriate colors
    ## according to 'name', and then use those to create the
    ## theme. The first 16 colors correspond to trellis.settings
    ## colors, the 17th is the background color.

    if (color)
    {
        ## color colors
        can.col <-
            if (name == "windows" || name == "X11")
                c("#000000", "#00ffff", "#ff00ff", "#00ff00",
                  "#ff7f00", "#007eff", "#ffff00", "#ff0000",
                  "#c6ffff", "#ffc3ff", "#c8ffc8", "#ffd18f",
                  "#a9e2ff", "#ffffc3", "#ff8c8a", "#aaaaaa",
                  "#909090")
            else if (name %in% c("postscript", "pdf", "xfig"))
                c("#000000", "#00ffff", "#ff00ff", "#00ff00",
                  "#ff7f00", "#0080ff", "#ffff00", "#ff0000",
                  "#ccffff", "#ffccff", "#ccffcc", "#ffe5cc",
                  "#cce6ff", "#ffffcc", "#ffcccc", "#e6e6e6",
                  "transparent")
            else ## default, same as X11 for now
                c("#000000", "#00FFFF", "#FF00FF", "#00FF00",
                  "#FF7F00", "#007EFF", "#FFFF00", "#FF0000",
                  "#C6FFFF", "#FFC3FF", "#C8FFC8", "#FFD18F",
                  "#A9E2FF", "#FFFFC3", "#FF8C8A", "#AAAAAA",
                  "#909090")
    }
    else ## b&w colors, same for all devices (8:16 actually unnecessary)
        can.col <-
            c("#000000", "#999999", "#4C4C4C", "#E6E6E6", "#F2F2F2",
              "#B2B2B2", "#000000", "#030303", "#050505", "#080808",
              "#0A0A0A", "#0D0D0D", "#0F0F0F", "#121212", "#151515",
              "#171717", "transparent")

    ## color settings, modified later if postscript or color = FALSE
    ans <-
        list(background = list(col = can.col[17]), 
             add.line = list(col = can.col[1], lty = 1, lwd = 1),
             add.text = list(cex = 1, col = can.col[1], font = 1),
             bar.fill = list(col = can.col[2]),
             box.dot = list(col = can.col[1], cex = 1, font = 1, pch =
             16),
             box.rectangle = list(col = can.col[2], lty = 1, lwd = 1),
             box.umbrella = list(col = can.col[2], lty = 2, lwd = 1),
             dot.line = list(col = can.col[16], lty = 1, lwd = 1),
             dot.symbol = list(cex = 0.8, col = can.col[2], font = 1,
             pch = 16),
             plot.line = list(col = can.col[2], lty = 1, lwd = 1),
             plot.symbol = list(cex = 0.8, col = can.col[2], font = 1,
             pch = 1),
             reference.line = list(col = can.col[16], lty = 1, lwd =
             1),
             strip.background = list(col = can.col[c(12, 11, 9, 13,
                                     10, 15, 14)]),
             strip.shingle = list(col = can.col[c(5, 4, 2, 6, 3, 8,
             7)]), superpose.line = list(col = can.col[2:8], lty =
             c(1, 1, 1, 1, 1, 1, 1), lwd = c(1, 1, 1, 1, 1, 1, 1)),
             regions = list(col = rev(cm.colors(100))),
             superpose.symbol = list(cex = c(0.8, 0.8, 0.8, 0.8, 0.8,
             0.8, 0.8), col = can.col[2:8], font = c(1, 1, 1, 1, 1, 1,
             1), pch = c("o", "o", "o", "o", "o", "o", "o")),
             axis.line = list(line = 0, col = can.col[1], lty = 1, lwd
             = 1),
             box.3d = list(col = can.col[1], lty = 1, lwd = 1),
             par.xlab.text = list(cex = 1, col = can.col[1], font =
             1),
             par.ylab.text = list(cex = 1, col = can.col[1], font =
             1),
             par.main.text = list(cex = 1.2, col = can.col[1], font =
             2),
             par.sub.text = list(cex = 1, col = can.col[1], font = 2))

    if (color) {
        if (name == "postscript" || name == "pdf") {
            ans$plot.symbol$col <- can.col[6]
            ans$plot.line$col <- can.col[6]
            ans$dot.symbol$col <- can.col[6]
            ans$box.rectangle$col <- can.col[6]
            ans$box.umbrella$col <- can.col[6]
            ans$superpose.symbol$col <- c(can.col[c(6, 3, 4, 8)],
                                          "orange", "darkgreen", "brown")
            ans$superpose.line$col <- c(can.col[c(6, 3, 4, 8)],
                                          "orange", "darkgreen", "brown")
        }
    }
    else {
        ## black and white settings
        ans$bar.fill$col <- can.col[5]
        ans$box.dot$col <- can.col[1]
        ans$box.rectangle$col <- can.col[1]
        ans$box.umbrella$col <- can.col[1]
        ans$box.umbrella$lty <- 2
        ans$dot.line$col <- can.col[4]
        ans$dot.symbol$col <- can.col[1]
        ans$dot.symbol$cex <- 0.85
        ans$plot.line$col <- can.col[1]
        ans$plot.symbol$col <- can.col[1]
        ans$regions$col <- gray(29:128/128)
        ans$reference.line$col <- can.col[4]
        ans$strip.background$col <- can.col[rep(5, 7)]
        ans$strip.shingle$col <- can.col[rep(6, 7)]
        ans$superpose.line$col <- can.col[rep(1, 7)]
        ans$superpose.line$lty <- 1:7
        ans$superpose.symbol$col <- can.col[rep(1, 7)]
        ans$superpose.symbol$cex <- rep(0.7, 7)
        ans$superpose.symbol$pch <- c(1,3,6,0,5,16,17)
        ##ans$superpose.symbol$pch <- c("o","+",">","s","w","#","{")
    }
    ans
}   





trellis.par.get <-
    function(name = NULL)
{
    ## this first block is needed to handle the unfortunate but
    ## inevitable case where a high level lattice function tries to
    ## access trellis parameters when in fact no device is open. There
    ## is of course no way of knowing at this point which device the
    ## final plot will be drawn to, so the default that would have
    ## been created using trtellis.device() without arguments is used

    if (!(.Device %in% names(lattice.theme))) {
        lattice.theme[[.Device]] <<- list()
        device <-
            if (is.null(dev.list()))
                getOption("device")
            else .Device
        color <- !(device == "postscript")
        lset(canonical.theme(device, color))
    }
    if (is.null(name))
        lattice.theme[[.Device]]
    else if (name %in% names(lattice.theme[[.Device]]))
        lattice.theme[[.Device]][[name]]
    else NULL
}

trellis.par.set <-
    function(name, value)
{
    ## if (name %in% names(lattice.theme[[.Device]])) NEEDED as a safeguard ?
    lattice.theme[[.Device]][[name]] <<- value
}



trellis.device <-
    function(device = getOption("device"),
             color = !(dev.name == "postscript"),
             theme = NULL,
             bg = NULL,
             new = TRUE,
             retain = FALSE,
             ...)
{
    ## Get device function
    if (is.character(device))
    {
        device.call <- get(device)
        dev.name <- device
    }
    else {
        device.call <- device
        dev.name <- deparse(substitute(device))
    }
    ## Start the new device if necessary.
    ## new = FALSE ignored if no devices open.
    if (new || .Device == "null device")
    {
        device.call(...)
        .lattice.print.more <<- FALSE
    }
    ## If retain = TRUE, retain existing settings for this
    ## particular device, if any. Ignores new.
    if (retain && .Device %in% names(lattice.theme))
    {
        if (!is.null(theme)) lset(theme) 
        return(invisible())
    }
    lset(canonical.theme(name=.Device, color=color))
    if (!is.null(theme)) lset(theme)
    if (!is.null(bg)) trellis.par.set("background", list(col = bg))
}




lset <- function(theme = col.whitebg())
{
    for (item in names(theme)) {
        foo <- trellis.par.get(item)
        bar <- theme[[item]]
        foo[names(bar)] <- bar
        trellis.par.set(item, foo)
    }
}







show.settings <- function(x = NULL)
{
    if (is.null(dev.list())) trellis.device()
    theme <- trellis.par.get()
    if (is.null(theme) && is.null(x)) print("No active device") ## shouldn't happen
    else {
        if (is.null(theme)) { ## also shouldn't happen
            print("No device is currently active but a theme has been explicitly specified.")
            print("Default device options will be used to fill missing components")
            print("of specified theme, if any.")
            theme <- canonical.theme(getOption("device"))
        }
        if (!is.null(x)) {
            for (item in names(x)) {
                foo <- x[[item]]
                theme[[item]][names(foo)] <- foo
            }
        }
        n.row <- 13
        n.col <- 9
        heights.x <- rep(1, n.row)
        heights.units <- rep("lines", n.row)
        heights.units[c(2, 5, 8, 11)] <- "null"
        widths.x <- rep(1, n.row)
        widths.units <- rep("lines", n.row)
        widths.units[c(2, 4, 6, 8)] <- "null"
        page.layout <- grid.layout(nrow = n.row, ncol = n.col,
                                   widths = unit(widths.x, widths.units),
                                   heights = unit(heights.x, heights.units))
        if (!.lattice.print.more) grid.newpage()
        .lattice.print.more <<- FALSE
        grid.rect(gp = gpar(fill = theme$background$col,
                  col = "transparent"))
        push.viewport(viewport(layout = page.layout))
        superpose.symbol <- theme$superpose.symbol
        len <- length(superpose.symbol$col)
        push.viewport(viewport(layout.pos.row = 2,
                               layout.pos.col = 2,
                               yscale = c(0,len+1),
                               xscale = c(0,len+1)))
        for (i in 1:len) {
            lpoints(y = rep(i, len), x = 1:len,
                    col = superpose.symbol$col[i],
                    font = superpose.symbol$font[i],
                    cex = superpose.symbol$cex[i],
                    pch = superpose.symbol$pch[i])
        }
        pop.viewport()
        grid.text(lab = "superpose.symbol",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 3, layout.pos.col = 2))
        superpose.line <- theme$superpose.line
        len <- length(superpose.line$col)
        push.viewport(viewport(layout.pos.row = 2,
                               layout.pos.col = 4,
                               yscale = c(0,len+1),
                               xscale = c(0,1)))
        for (i in 1:len) {
            llines(y = rep(i, 2), x = c(0,1),
                   col = superpose.line$col[i],
                   lty = superpose.line$lty[i],
                   lwd = superpose.line$lwd[i])
        }
        pop.viewport()
        grid.text(lab = "superpose.line",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 3, layout.pos.col = 4))
        strip.background <- theme$strip.background
        len <- length(strip.background$col)
        push.viewport(viewport(layout.pos.row = 2,
                               layout.pos.col = 6,
                               yscale = c(0,len+1),
                               xscale = c(0,1)))
        for (i in 1:len) {
            grid.rect(y = unit(i, "native"), h = unit(.5, "native"),
                      gp = gpar(fill = strip.background$col[i]))
        }
        pop.viewport()
        grid.text(lab = "strip.background",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 3, layout.pos.col = 6))
        strip.shingle <- theme$strip.shingle
        len <- length(strip.shingle$col)
        push.viewport(viewport(layout.pos.row = 2,
                               layout.pos.col = 8,
                               yscale = c(0,len+1),
                               xscale = c(0,1)))
        for (i in 1:len) {
            grid.rect(y = unit(i, "native"), h = unit(.5, "native"),
                      gp = gpar(fill = strip.shingle$col[i]))
        }
        pop.viewport()
        grid.text(lab = "strip.shingle",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 3, layout.pos.col = 8))
        push.viewport(viewport(layout.pos.row = 5,
                               layout.pos.col = 2,
                               yscale = extend.limits(c(0,6)),
                               xscale = c(0,6)))
        x <- c(1,2,3,4,5)
        dot.line <- theme$dot.line
        dot.symbol <- theme$dot.symbol
        panel.abline(h=1:5, col=dot.line$col,
                     lty=dot.line$lty, lwd=dot.line$lwd)
        panel.xyplot(x = x, y = x, col = dot.symbol$col, pch = dot.symbol$pch)
        grid.rect()
        pop.viewport()
        grid.text(lab = "dot.[symbol, line]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 6, layout.pos.col = 2))
        box.rectangle <- theme$box.rectangle
        box.dot <- theme$box.dot
        box.umbrella <- theme$box.umbrella
        push.viewport(viewport(layout.pos.row = 5,
                               layout.pos.col = 4,
                               yscale = c(-1.5,1.5),
                               xscale = c(0,6)))
        push.viewport(viewport(height = unit(.1, "npc")))
        grid.rect(width = 1/3, 
                  gp = gpar(lwd = box.rectangle$lwd, 
                  lty = box.rectangle$lty, col = box.rectangle$col))
        grid.lines(x = unit(c(1/6, 1/3), "npc"), 
                   y = unit(c(0.5, 0.5), "npc"),
                   gp = gpar(col = box.umbrella$col,
                   lwd = box.umbrella$lwd, lty = box.umbrella$lty))
        grid.lines(x = unit(c(2/3, 5/6), "npc"), 
                   y = unit(c(0.5, 0.5), "npc"),
                   gp = gpar(col = box.umbrella$col,
                   lwd = box.umbrella$lwd, lty = box.umbrella$lty))
        grid.lines(x = unit(rep(1/6, 2), "npc"), 
                   y = unit(c(0, 1), "npc"),
                   gp = gpar(col = box.umbrella$col, 
                   lwd = box.umbrella$lwd, lty = box.umbrella$lty))
        grid.lines(x = unit(rep(5/6, 2), "npc"), 
                   y = unit(c(0, 1), "npc"),
                   gp = gpar(col = box.umbrella$col, 
                   lwd = box.umbrella$lwd, lty = box.umbrella$lty))
        if (is.character(box.dot$pch)) 
            grid.text(x = unit(0.5, "npc"), y = unit(0.5, "npc"),
                      lab = box.dot$pch,
                      gp = gpar(col = box.dot$col, cex = box.dot$cex))
        else
            grid.points(x = unit(0.5, "npc"), y = unit(0.5, "npc"),
                        pch = box.dot$pch, size = unit(box.dot$cex * 2.5, "mm"),
                        gp = gpar(col = box.dot$col))
        pop.viewport()
        grid.rect()
        pop.viewport()
        grid.text(lab = "box.[dot, rectangle, umbrella]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 6, layout.pos.col = 4))
        add.text <- theme$add.text
        add.line <- theme$add.line
        push.viewport(viewport(layout.pos.row = 5,
                               layout.pos.col = 6,
                               yscale = c(-1,1),
                               xscale = c(0,1)))
        x <- seq(.1, .9, length = 50)
        y <- .9 * sin(.1+11*x)
        llines(x = x, y = y, type = "l", col = add.line$col,
               lty = add.line$lty, lwd = add.line$lwd)
        ltext(lab = c("Hello", "World"),
              x = c(.25, .75), y = c(-.5, .5),
              col = add.text$col, cex = add.text$cex, font = add.text$font)
        grid.rect()
        pop.viewport()
        grid.text(lab = "add.[line, text]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 6, layout.pos.col = 6))
        reference.line <- theme$reference.line
        push.viewport(viewport(layout.pos.row = 5,
                               layout.pos.col = 8,
                               yscale = c(0,4),
                               xscale = c(0,4)))
        panel.grid(col = reference.line$col,
                   lwd = reference.line$lwd,
                   lty = reference.line$lty)
        grid.rect()
        pop.viewport()
        grid.text(lab = "reference.line",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 6, layout.pos.col = 8))
        plot.symbol <- theme$plot.symbol
        plot.line <- theme$plot.line
        push.viewport(viewport(layout.pos.row = 8,
                               layout.pos.col = 2,
                               yscale = c(-1.1,1.1),
                               xscale = c(-.1,1.1)))
        x <- seq(.1, .9, length = 20)
        y <- .9 * sin(.1+11*x)
        panel.xyplot(x = x+.05, y = y+.1, type = "l",
                     lty = plot.line$lty,
                     col = plot.line$col,
                     lwd = plot.line$lwd)
        panel.xyplot(x = x-.05, y = y-.1,
                     col = plot.symbol$col,
                     font = plot.symbol$font,
                     pch = plot.symbol$pch,
                     cex = plot.symbol$cex)
        grid.rect()
        pop.viewport()
        grid.text(lab = "plot.[symbol, line]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 9, layout.pos.col = 2))
        bar.fill <- theme$bar.fill
        push.viewport(viewport(layout.pos.row = 8,
                               layout.pos.col = 4,
                               yscale = extend.limits(c(0,6)),
                               xscale = extend.limits(c(1,10))))
        grid.rect(x = c(3.5, 4.5, 5.5, 6.5, 7.5), w = rep(5,5),
                  y = c(1,2,3,4,5), h = rep(.5, ,5),
                  default.units = "native",
                  gp = gpar(fill = bar.fill$col))
        grid.rect()
        pop.viewport()
        grid.text(lab = "plot.shingle[bar.fill]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 9, layout.pos.col = 4))
        push.viewport(viewport(layout.pos.row = 8,
                               layout.pos.col = 6,
                               yscale = extend.limits(c(0,7)),
                               xscale = extend.limits(c(0,7))))
        grid.rect(y = c(.5, 1, 1.5, 2, 2.5, 3, 3.5), w = rep(1,7),
                  x = c(.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5), h = 1:7,
                  default.units = "native",
                  gp = gpar(fill = bar.fill$col))
        grid.rect()
        pop.viewport()
        grid.text(lab = "histogram[bar.fill]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 9, layout.pos.col = 6))
        push.viewport(viewport(layout.pos.row = 8,
                               layout.pos.col = 8,
                               yscale = extend.limits(c(0,6)),
                               xscale = c(0,7)))
        grid.rect(x = rev(c(.5, 1, 1.5, 2, 2.5, 3)), h = rep(.5, 6),
                  y = c(.5, 1.5, 2.5, 3.5, 4.5, 5.5), w = 6:1,
                  default.units = "native",
                  gp = gpar(fill = bar.fill$col))
        grid.rect()
        pop.viewport()
        grid.text(lab = "barchart[bar.fill]",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 9, layout.pos.col = 8))
        regions <- theme$regions
        len <- length(regions$col)
        push.viewport(viewport(layout.pos.row = 11,
                               layout.pos.col = 2,
                               xscale = c(0,len+1)))
        for (i in 1:len)
            grid.rect(x = i, w = 1, default.units = "native",
                      gp = gpar(col = NULL,  fill = regions$col[i]))
        grid.rect()
        pop.viewport()
        grid.text(lab = "regions",
                  gp = gpar(fontsize = 8),
                  vp = viewport(layout.pos.row = 12, layout.pos.col = 2))
    }    
    invisible()
}










    


### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA





prepanel.default.splom <-
    function(x, y, type, ...)
{
    list(xlim = c(0,1),
         ylim = c(0,1),
         dx = 1,
         dy = 1)
}




panel.splom <-
    function(...)
    panel.xyplot(...)





panel.pairs <-
    function(z, panel = "panel.splom", groups = NULL,
             panel.subscripts,
             subscripts,
             fontsize.small = 8,
             pscales = 5,
             ...)
{
    panel <- 
        if (is.function(panel)) panel 
        else if (is.character(panel)) get(panel)
        else eval(panel)

    axis.line <- trellis.par.get("axis.line")
    n.var <- ncol(z)

    if(n.var>0) {
        ## there must be a better way to do the foll:
        lim <- list(1:n.var)
        for(i in 1:n.var) {
            lim[[i]] <- extend.limits(range(as.numeric(z[,i])))
        }
        ## should be further complicated by allowing for customization by
        ## prepanel functions --- prepanel(z[i], z[j]) etc
    }
    ## maybe (ideally) this should be affected by scales

    draw <- is.numeric(pscales) && pscales!=0 # whether axes to be drawn

    splom.layout <- grid.layout(nrow=n.var, ncol=n.var)

    if(n.var > 0 && any(subscripts)) {

        push.viewport(viewport(layout=splom.layout))

        for(i in 1:n.var)
            for(j in 1:n.var)
            {
                push.viewport(viewport(layout.pos.row = n.var-i+1,
                                       layout.pos.col = j,
                                       clip = TRUE,
                                       gp = gpar(fontsize = fontsize.small),
                                       xscale = lim[[j]],
                                       yscale = lim[[i]]))

                if(i == j)
                {
                    if (!is.null(colnames(z)))
                        grid.text(colnames(z)[i],
                                  gp = gpar(fontsize = 10))
                    if (draw) {
                        ## plot axes

                        if (is.factor(z[,i])) {
                            axls <- 1:nlevels(z[,i])
                            nal <- length(axls)/2+.5

                            for(tt in seq(along=axls)) {
                                if(tt <= nal) {
                                    
                                    grid.lines(y = unit(rep(axls[tt],2), "native"),
                                               x = unit(c(1,1),"npc") - unit(c(0,.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = levels(z[,i])[tt],
                                              x = unit(1,"npc") - unit(.5, "lines"),
                                              y = unit(axls[tt], "native"),
                                              just = c("right", "centre"))
                                    
                                    grid.lines(x = unit(rep(axls[tt],2), "native"),
                                               y = unit(c(0,.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = levels(z[,i])[tt],
                                              rot = 90,
                                              y = unit(0.5, "lines"),
                                              x = unit(axls[tt], "native"),
                                              just = c("bottom", "centre"))
                                    
                                }
                                if(tt >=nal) {
                                    
                                    grid.lines(y = unit(rep(axls[tt],2), "native"),
                                               x = unit(c(0,0.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = levels(z[,i])[tt],
                                              x = unit(0.5, "lines"),
                                              y = unit(axls[tt], "native"),
                                              just = c("left", "centre"))
                                    
                                    grid.lines(x = unit(rep(axls[tt],2), "native"),
                                               y = unit(c(1,1),"npc") - unit(c(0,.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = levels(z[,i])[tt], rot = 90,
                                              y = unit(1,"npc") - unit(.5, "lines"),
                                              x = unit(axls[tt], "native"),
                                              just = c("top", "centre"))
                                    
                                }
                                
                            }
                            
                        }
                        else {
                        
                            axls <-
                                if (is.list(pscales) && !is.null(pscales[[i]]$at))
                                    pscales[[i]]$at
                                else
                                    lpretty(lim[[i]], n = pscales)

                            labels <-
                                if (is.list(pscales) && !is.null(pscales[[i]]$lab))
                                    pscales[[i]]$lab
                            ## should be rendered like factors ?
                                else
                                    as.character(axls)

                            axid <- axls>lim[[i]][1] & axls <lim[[i]][2]
                            axls <- axls[axid]
                            labels <- labels[axid]
                            nal <- length(axls)/2+.5

                            for(tt in seq(along=axls)) {
                                if(tt <= nal) {
                                    
                                    grid.lines(y = unit(rep(axls[tt],2), "native"),
                                               x = unit(c(1,1),"npc") - unit(c(0,.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = as.character(axls[tt]),
                                              x = unit(1,"npc") - unit(.5, "lines"),
                                              y = unit(axls[tt], "native"),
                                              just = c("right", "centre"))
                                    
                                    grid.lines(x = unit(rep(axls[tt],2), "native"),
                                               y = unit(c(0,.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = as.character(axls[tt]),
                                              y = unit(0.5, "lines"),
                                              x = unit(axls[tt], "native"),
                                              just = c("centre", "left"))
                                    
                                }
                                if(tt >=nal) {
                                    
                                    grid.lines(y = unit(rep(axls[tt],2), "native"),
                                               x = unit(c(0,0.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = as.character(axls[tt]),
                                              x = unit(0.5, "lines"),
                                              y = unit(axls[tt], "native"),
                                              just = c("left", "centre"))
                                    
                                    grid.lines(x = unit(rep(axls[tt],2), "native"),
                                               y = unit(c(1,1),"npc") - unit(c(0,.25), "lines"),
                                               gp = gpar(col = axis.line$col))
                                    
                                    grid.text(label = as.character(axls[tt]),
                                              y = unit(1,"npc") - unit(.5, "lines"),
                                              x = unit(axls[tt], "native"),
                                              just = c("centre", "right"))
                                    
                                }
                                
                            }
                        }    
                    }

                    grid.rect()

                }
                else
                {
                    if(!panel.subscripts)
                        panel(x=as.numeric(z[subscripts, j]),
                              y=as.numeric(z[subscripts, i]), ...)

                    else panel(x=as.numeric(z[subscripts, j]),
                               y=as.numeric(z[subscripts, i]),
                               groups = groups,
                               subscripts = subscripts, ...)

                    grid.rect()
                }
                pop.viewport()
            }
        pop.viewport()
    }
}




splom <-
    function(formula,
             data = parent.frame(),
             aspect = 1,
             between = list(x = 0.5, y = 0.5),
             layout = NULL,
             panel = "panel.splom",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab = "Scatter Plot Matrix",
             xlim,
             ylab = NULL,
             ylim,
             superpanel = "panel.pairs",
             pscales = 5,
             varnames,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ## dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing
    


    formname <- deparse(substitute(formula))
    formula <- eval(substitute(formula), data, parent.frame())

    form <-
        if (inherits(formula, "formula"))
            latticeParseFormula(formula, data)
        else 
            list(left = NULL,
                 right = as.data.frame(formula),
                 condition = NULL,
                 left.name = "",
                 right.name = formname)


    ##form <- latticeParseFormula(formula, data)

    cond <- form$condition


    number.of.cond <- length(cond)
    x <- as.data.frame(form$right)
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, nrow(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }
    if (!missing(varnames)) colnames(x) <-
        eval(substitute(varnames), data, parent.frame())
    
    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    subscr <- seq(along=x[,1])
    x <- x[subset,, drop = TRUE]
    subscr <- subscr[subset, drop = TRUE]
    
    ##if(!(is.numeric(x) && is.numeric(y)))
    ##    warning("Both x and y should be numeric")    WHAT ?


    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          between = between,
                          panel = superpanel,
                          strip = strip,
                          xlab = xlab,
                          ylab = ylab), dots))

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- "Scatter Plot Matrix"
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab <- NULL

    ## Step 2: Compute scales.common (leaving out limits for now)

    ## It is not very clear exactly what effect scales is supposed
    ## to have. Not much in Trellis (probably), but there are certain
    ## components which are definitely relevant, and certail others
    ## (like log) which can be used in innovative ways. However, I'm
    ## postponing all that to later, if at all,and for now TOTALLY
    ## ignoring scales
    
    ##scales <- eval(substitute(scales), data, parent.frame())
    ##if (is.character(scales)) scales <- list(relation = scales)
    scales <- list(relation = "same", draw = FALSE)
    foo <- c(foo, 
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }
    if (have.xlim || have.ylim) {
        warning("Limits cannot be explicitly specified")
    }
    have.xlim <- TRUE
    have.ylim <- TRUE
    xlim <- c(0,1)
    ylim <- c(0,1)
    
    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
#    if (have.xlog) {
#        xlog <- foo$x.scales$log
#        xbase <-
#            if (is.logical(xlog)) 10
#            else if (is.numeric(xlog)) xlog
#            else if (xlog == "e") exp(1)
#
#        x <- log(x, xbase)
#        if (have.xlim) xlim <- log(xlim, xbase)
#    }
#    if (have.ylog) {
#        ylog <- foo$y.scales$log
#        ybase <-
#            if (is.logical(ylog)) 10
#            else if (is.numeric(ylog)) ylog
#            else if (ylog == "e") exp(1)
#
#        y <- log(y, ybase)
#        if (have.ylim) ylim <- log(ylim, ybase)
#    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- FALSE
    for (j in 1:ncol(x))
        id.na <- id.na | is.na(x[,j])
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args


    foo$panel.args.common <-
        c(list(z = x,
               panel = panel,
               panel.subscripts = subscripts,
               groups = groups, # xscales = foo$x.scales, yscales = foo$y.scales,
               pscales = pscales),
          dots)

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >=
                              levels(var)[[cond.current.level[i]]][1])
                             & (var <=
                                levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }

                    foo$panel.args[[panel.number]] <-
                        list(subscripts = subscr[id])

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.splom,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}



### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA



prepanel.default.tmd <-
    function(...)
    prepanel.default.xyplot(...)



panel.tmd <- function(...) {
    panel.abline(h=0)
    panel.xyplot(...)
}


## Fixme: log scales not handled
tmd <-
    function(object,
             aspect = "fill",
             as.table = object$as.table,
             between = list(x=object$x.between,y=object$y.between),
             key = object$key,
             layout = object$layout,
             main = object$main,
             page = object$page,
             panel = "panel.tmd",
             par.strip.text = object$par.strip.text, 
             prepanel = NULL,
             scales = list(),
             strip = object$strip,
             sub = object$sub,
             xlab = "mean",
             xlim = NULL,
             ylab = "difference",
             ylim = NULL,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(as.table = as.table,
                          aspect = aspect,
                          between = between,
                          key = key,
                          page = page,
                          main = main,
                          panel = panel,
                          sub = sub,
                          par.strip.text = par.strip.text,
                          strip = strip,
                          xlab = xlab,
                          ylab = ylab), dots))
                          

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- form$right.name
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <- form$left.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo, 
             do.call("construct.scales", scales))

    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) { ## problem
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) { ## problem
        ylog <- foo$y.scales$log
        ybase <-
            if (is.logical(ylog)) 10
            else if (is.numeric(ylog)) ylog
            else if (ylog == "e") exp(1)

        y <- log(y, ybase)
        if (have.ylim) ylim <- log(ylim, ybase)
    }
    
    ## Step 5: Process cond

    foo$condlevels <- object$condlevels

    ## Step 6: Evaluate layout, panel.args.common and panel.args

    foo$panel.args.common <- c(object$panel.args.common, dots)

    if (!missing(layout)) {
        number.of.cond <- length(foo$condlevels)
        cond.max.level <- integer(number.of.cond)
        for(i in 1:number.of.cond) {
            cond.max.level[i] <-
                if (is.character(foo$condlev[[i]])) length(foo$condlev[[i]])
                else nrow(foo$condlev[[i]])
        }
        foo$skip <- !unlist(lapply(object$panel.args , is.list))
        layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    }
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- object$panel.args

    if ("x" %in% names(foo$panel.args.common)) {
        ## this would happen with subscripts. assuming that
        ## y would also be there then
        q < foo$panel.args.common
        x <- (q$x+q$y)/2
        y <- q$y-q$x       # will stop if any errors, not putting any more handlers
        foo$panel.args.common$x <- x
        foo$panel.args.common$y <- y
    }
    else {
        count <- 1
        for (p in foo$panel.args)
            if (is.logical(p)) # which means skip = T for this panel
                count <- count + 1 
            else {
                x <- (p$x+p$y)/2
                y <- p$y-p$x

                foo$panel.args[[count]]$x <- x
                foo$panel.args[[count]]$y <- y

                count <- count + 1
            }
    }    

    foo <- c(foo,
             limits.and.aspect(prepanel.default.tmd,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}


update.trellis <-
    function(object,
             aspect,
             as.table,
             between,
             key,
             layout,
             main,
             page,
             panel,
             par.strip.text,
             scales,
             skip,
             strip,
             sub,
             xlab,
             xlim,
             ylab,
             ylim,
             ...)
{
    
    if(!missing(aspect)){
        if (is.numeric(aspect)){
            object$aspect.ratio <- aspect
            object$aspect.fill <- FALSE
        }
        else if (is.character(aspect) && aspect=="fill")
            object$aspect.fill <- TRUE
        else warning("Inappropriate value of aspect")
    }
    
    if(!missing(as.table)){
        if (is.logical(as.table)) object$as.table <- as.table
        else warning("Inappropriate value of as.table")
    }
    
    if(!missing(between)){
        if ("x" %in% names(between)) object$x.between <- between$x
        if ("y" %in% names(between)) object$y.between <- between$y
    }
    
    if(!missing(key)){
        object$key <- key
    }
    
    if(!missing(layout)){
        if (length(layout)==2){
            object$layout[3] <- ceiling(object$layout[1]*object$layout[2]*object$layout[3]/
                                     max(layout[1]*layout[2], layout[2]))
            object$layout[1] <- layout[1]
            object$layout[2] <- layout[2]
        }
        else if (length(layout)==3){
            object$layout <- layout
        }
        else warning("Inappropriate value of layout")
    }



    
    if(!missing(main)){
        if (is.character(main))
            if (is.null(object$main)) object$main <-
                list(label = main, col = "black", cex = 1, font = 1)
            else object$main$label <- main

        else if (is.list(main)) {
            object$main[names(main)] <- main
        }
        else if (is.null(main)) object$main <- NULL
    }
    
    if(!missing(sub)){
        if (is.character(sub))
            if (is.null(object$sub)) object$sub <-
                list(label = sub, col = "black", cex = 1, font = 1)
            else object$sub$label <- sub

        else if (is.list(sub)) {
            object$sub[names(sub)] <- sub
        }
        else if (is.null(sub)) object$sub <- NULL
    }
    
    if(!missing(xlab)){
        if (is.character(xlab))
            if (is.null(object$xlab)) object$xlab <-
                list(label = xlab, col = "black", cex = 1, font = 1)
            else object$xlab$label <- xlab

        else if (is.list(xlab)) {
            object$xlab[names(xlab)] <- xlab
        }
        else if (is.null(xlab)) object$xlab <- NULL
    }
    
    if(!missing(ylab)){
        if (is.character(ylab))
            if (is.null(object$ylab)) object$ylab <-
                list(label = ylab, col = "black", cex = 1, font = 1)
            else object$ylab$label <- ylab

        else if (is.list(ylab)) {
            object$ylab[names(ylab)] <- ylab
        }
        else if (is.null(ylab)) object$ylab <- NULL
    }
    

    
    if(!missing(page)){
        object$page <- page
    }
    
    if(!missing(panel)){
        panel <- 
            if (is.function(panel)) panel 
            else if (is.character(panel)) get(panel)
            else eval(panel)

        if (object$fname == "splom")
            object$panel.args.common$panel <- panel
        else object$panel <- panel
    }
    
    if(!missing(par.strip.text)){
        if (is.list(par.strip.text))
            object$par.strip.text[names(par.strip.text)] <-
                par.strip.text
        else warning("par.strip.text must be a list")
    }
    
    if(!missing(skip)){
        warning("sorry, but skip cannot be changed by update")
    }
    
    if(!missing(strip)){
        if (is.logical(strip)){
            if (strip) object$strip <- strip.default
        else object$strip <- FALSE
        }
        else object$strip <- strip
    }
    
    if(!missing(xlim)){
        if (!is.list(object$x.limits)) object$x.limits <- xlim
        else warning("xlim cannot be specified unless relation = same")
    }
    
    if(!missing(ylim)){
        if (!is.list(object$y.limits)) object$y.limits <- ylim
        else warning("ylim cannot be specified unless relation = same")
    }

    if(!missing(scales)){
        
        if ("relation" %in% names(scales))
            warning("relation cannot be changed via update")
        
      

      
        if ("alternating" %in% names(scales)){
            if (is.logical(scales$alternating))
                if (scales$alternating){
                    object$x.scales$alternating <- c(1,2)
                    object$y.scales$alternating <- c(1,2)
                }
                else {
                    object$x.scales$alternating <- 1
                    object$y.scales$alternating <- 1
                }
            else if (is.numeric(scales$alternating)) {
                object$x.scales$alternating <- scales$alternating
                object$y.scales$alternating <- scales$alternating
            }
        }
        
        
        
        
        if ("x" %in% names(scales) &&
            "alternating" %in% names(scales$x)){
            if (is.logical(scales$x$alternating))
                if (scales$x$alternating){
                    object$x.scales$alternating <- c(1,2)
                }
                else {
                    object$x.scales$alternating <- 1
                }
            else if (is.numeric(scales$x$alternating)) {
                object$x.scales$alternating <- scales$x$alternating
            }
        }
        
        
        
        
        
        if ("y" %in% names(scales) &&
            "alternating" %in% names(scales$y)){
            if (is.logical(scales$y$alternating))
                if (scales$y$alternating){
                    object$y.scales$alternating <- c(1,2)
                }
                else {
                    object$y.scales$alternating <- 1
                }
            else if (is.numeric(scales$y$alternating)) {
                object$y.scales$alternating <- scales$y$alternating
            }
        }


        object$x.scales[names(scales)] <- scales
        if ("x" %in% names(scales)){
                
            object$x.scales[names(scales$x)] <- scales$x
                
            if ("relation" %in% names(scales$x))
                warning("relation cannot be changed via update")
        }
        
        object$y.scales[names(scales)] <- scales
        if ("y" %in% names(scales)){
                
            object$y.scales[names(scales$y)] <- scales$y
                
            if ("relation" %in% names(scales$y))
                warning("relation cannot be changed via update")
        }

    }

    object
}


### Copyright 2001  Deepayan Sarkar <deepayan@stat.wisc.edu>
###
### This file is part of the lattice library for R.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA




prepanel.default.xyplot <-
    function(x, y, type, ...)
{
    ord <- order(x)
    list(xlim = range(x),
         ylim = range(y),
         dx = diff(x[ord]),
         dy = diff(y[ord]))
}




panel.xyplot <-
    function(x, y, type="p",
             pch = plot.symbol$pch,
             col,
             col.line = plot.line$col,
             col.symbol = plot.symbol$col,
             lty = plot.line$lty,
             cex = plot.symbol$cex,
             lwd = plot.line$lwd, ...)
{
    if (length(x)>0) {


        if (!missing(col)) {
            if (missing(col.line)) col.line <- col
            if (missing(col.symbol)) col.symbol <- col
        }

        plot.symbol <- trellis.par.get("plot.symbol")
        plot.line <- trellis.par.get("plot.line")

        if ("o" %in% type || "b" %in% type)
            type <- c(type, "p", "l")


        if ("p" %in% type)
            lpoints(x = x, y = y, cex = cex,
                    col = col.symbol, pch=pch)


        if ("l" %in% type)
            llines(x=x, y=y, lty=lty, col=col.line, lwd=lwd)


        if ("h" %in% type)
            llines(x=x, y=y, type = "h", lty=lty, col=col.line, lwd=lwd)


        if ("s" %in% type) {
            ord <- order(x)
            n <- length(x)
            xx <- numeric(2*n-1)
            yy <- numeric(2*n-1)

            xx[2*1:n-1] <- x[ord]
            yy[2*1:n-1] <- y[ord]
            xx[2*1:(n-1)] <- x[ord][-1]
            yy[2*1:(n-1)] <- y[ord][-n]
            llines(x=xx, y=yy,
                   lty=lty, col=col.line, lwd=lwd)
        }
        if ("S" %in% type) {
            ord <- order(x)
            n <- length(x)
            xx <- numeric(2*n-1)
            yy <- numeric(2*n-1)

            xx[2*1:n-1] <- x[ord]
            yy[2*1:n-1] <- y[ord]
            xx[2*1:(n-1)] <- x[ord][-n]
            yy[2*1:(n-1)] <- y[ord][-1]
            llines(x=xx, y=yy,
                   lty=lty, col=col.line, lwd=lwd)
        }
        if ("r" %in% type) {
            panel.lmline(x, y, col = col.line, lty = lty, lwd = lwd, ...)
        }
        if ("smooth" %in% type) {
            panel.loess(x, y, col = col.line, lty = lty, lwd = lwd, ...)
        }
    }
}





xyplot <-
    function(formula,
             data = parent.frame(),
             aspect = "fill",
             layout = NULL,
             panel = "panel.xyplot",
             prepanel = NULL,
             scales = list(),
             strip = TRUE,
             groups = NULL,
             xlab,
             xlim,
             ylab,
             ylim,
             ...,
             subscripts = !is.null(groups),
             subset = TRUE)
{

    ##dots <- eval(substitute(list(...)), data, parent.frame())
    dots <- list(...)

    if (!is.function(panel)) panel <- eval(panel)
    if (!is.function(strip)) strip <- eval(strip)

    prepanel <-
        if (is.function(prepanel)) prepanel 
        else if (is.character(prepanel)) get(prepanel)
        else eval(prepanel)

    ## Step 1: Evaluate x, y, etc. and do some preprocessing

    form <- latticeParseFormula(formula, data)
    cond <- form$condition
    number.of.cond <- length(cond)
    y <- form$left
    x <- form$right
    if (number.of.cond == 0) {
        strip <- FALSE
        cond <- list(as.factor(rep(1, length(x))))
        layout <- c(1,1,1)
        number.of.cond <- 1
    }

    groups <- eval(substitute(groups), data, parent.frame())
    subset <- eval(substitute(subset), data, parent.frame())
    if ("subscripts" %in% names(formals(eval(panel)))) subscripts <- TRUE
    if(subscripts) subscr <- seq(along=x)
    x <- x[subset, drop = TRUE]
    y <- y[subset, drop = TRUE]
    if (subscripts) subscr <- subscr[subset, drop = TRUE]
    
    if (missing(xlab)) xlab <- form$right.name
    if (missing(ylab)) ylab <- form$left.name

    if(!(is.numeric(x) && is.numeric(y)))
        warning("Both x and y should be numeric")
    x <- as.numeric(x)
    y <- as.numeric(y)
    ## create a skeleton trellis object with the
    ## less complicated components:

    foo <- do.call("trellis.skeleton",
                   c(list(aspect = aspect,
                          strip = strip,
                          panel = panel,
                          xlab = xlab,
                          ylab = ylab), dots))
                          

    dots <- foo$dots # arguments not processed by trellis.skeleton
    foo <- foo$foo
    foo$call <- match.call()
    foo$fontsize.normal <- 10
    foo$fontsize.small <- 8

    ## This is for cases like xlab/ylab = list(cex=2)
    if (is.list(foo$xlab) && !is.character(foo$xlab$label))
        foo$xlab$label <- form$right.name
    if (is.list(foo$ylab) && !is.character(foo$ylab$label))
        foo$ylab$label <- form$left.name

    ## Step 2: Compute scales.common (leaving out limits for now)

    ##scales <- eval(substitute(scales), data, parent.frame())
    if (is.character(scales)) scales <- list(relation = scales)
    foo <- c(foo, 
             do.call("construct.scales", scales))


    ## Step 3: Decide if limits were specified in call:

    have.xlim <- !missing(xlim)
    if (!is.null(foo$x.scales$limit)) {
        have.xlim <- TRUE
        xlim <- foo$x.scales$limit
    }
    have.ylim <- !missing(ylim)
    if (!is.null(foo$y.scales$limit)) {
        have.ylim <- TRUE
        ylim <- foo$x.scales$limit
    }

    ## Step 4: Decide if log scales are being used:

    have.xlog <- !is.logical(foo$x.scales$log) || foo$x.scales$log
    have.ylog <- !is.logical(foo$y.scales$log) || foo$y.scales$log
    if (have.xlog) {
        xlog <- foo$x.scales$log
        xbase <-
            if (is.logical(xlog)) 10
            else if (is.numeric(xlog)) xlog
            else if (xlog == "e") exp(1)

        x <- log(x, xbase)
        if (have.xlim) xlim <- log(xlim, xbase)
    }
    if (have.ylog) {
        ylog <- foo$y.scales$log
        ybase <-
            if (is.logical(ylog)) 10
            else if (is.numeric(ylog)) ylog
            else if (ylog == "e") exp(1)

        y <- log(y, ybase)
        if (have.ylim) ylim <- log(ylim, ybase)
    }
    
    ## Step 5: Process cond

    cond <- lapply(cond, as.factorOrShingle, subset, drop = TRUE)
    cond.max.level <- unlist(lapply(cond, nlevels))


    id.na <- is.na(x)|is.na(y)
    for (var in cond)
        id.na <- id.na | is.na(var)
    if (!any(!id.na)) stop("nothing to draw")
    ## Nothing simpler ?

    foo$condlevels <- lapply(cond, levels)

    ## Step 6: Evaluate layout, panel.args.common and panel.args


    foo$panel.args.common <- dots
    if (subscripts) foo$panel.args.common$groups <- groups

    layout <- compute.layout(layout, cond.max.level, skip = foo$skip)
    plots.per.page <- max(layout[1] * layout[2], layout[2])
    number.of.pages <- layout[3]
    foo$skip <- rep(foo$skip, length = plots.per.page)
    foo$layout <- layout
    nplots <- plots.per.page * number.of.pages

    foo$panel.args <- as.list(1:nplots)
    cond.current.level <- rep(1,number.of.cond)
    panel.number <- 1 # this is a counter for panel number
    for (page.number in 1:number.of.pages)
        if (!any(cond.max.level-cond.current.level<0))
            for (plot in 1:plots.per.page) {

                if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
                else if(!any(cond.max.level-cond.current.level<0)) {

                    id <- !id.na
                    for(i in 1:number.of.cond)
                    {
                        var <- cond[[i]]
                        id <- id &
                        if (is.shingle(var))
                            ((var >= levels(var)[[cond.current.level[i]]][1])
                             & (var <= levels(var)[[cond.current.level[i]]][2]))
                        else (as.numeric(var) == cond.current.level[i])
                    }
                    foo$panel.args[[panel.number]] <-
                        list(x = x[id], y = y[id])
                    if (subscripts)
                        foo$panel.args[[panel.number]]$subscripts <-
                            subscr[id]

                    cond.current.level <-
                        cupdate(cond.current.level,
                                cond.max.level)
                }

                panel.number <- panel.number + 1
            }

    foo <- c(foo,
             limits.and.aspect(prepanel.default.xyplot,
                               prepanel = prepanel, 
                               have.xlim = have.xlim, xlim = xlim, 
                               have.ylim = have.ylim, ylim = ylim, 
                               x.relation = foo$x.scales$relation,
                               y.relation = foo$y.scales$relation,
                               panel.args.common = foo$panel.args.common,
                               panel.args = foo$panel.args,
                               aspect = aspect,
                               nplots = nplots))

    class(foo) <- "trellis"
    foo
}








### loading required library grid



## Need global variable to handle more in print.trellis
.lattice.print.more <- FALSE

autoload("loess.smooth", "modreg")
autoload("loess", "modreg")
library(grid)

.First.lib <- function(lib, pkg) {
  library.dynam( "lattice", pkg, lib )
}


