plsmo <-
  function(x, y, method=c("lowess", "supsmu", "raw"),
           xlab, ylab, add=FALSE, lty=1 : lc, col=par('col'),
           lwd=par('lwd'), iter=if(length(unique(y)) > 2) 3 else 0,
           bass=0, f=2 / 3, trim, fun, group=rep(1, length(x)),
           prefix, xlim, ylim, 
           label.curves=TRUE, datadensity=FALSE, scat1d.opts=NULL,
           lines.=TRUE, subset=TRUE, grid=FALSE, evaluate=NULL, ...)
{
  gfun <- ordGridFun(grid)
  nam <- as.character(sys.call())[2 : 3]
  method <- match.arg(method)
  if(missing(ylab))
    ylab <- label(y, units=TRUE, plot=TRUE, default=nam[2])
  Y <- as.matrix(y)
  p <- ncol(Y)
  if(!missing(subset)) {
    x     <- x[subset]
    Y     <- Y[subset,, drop=FALSE]
    group <- group[subset]
  }
    
  group <- as.factor(group)
  if(!missing(prefix)) levels(group) <- paste(prefix, levels(group))
  
  group <- as.factor(group)
  nna <- !(is.na(x) | (rowSums(is.na(Y)) == p) | is.na(group))
  x     <- x[nna]
  Y     <- Y[nna,, drop=FALSE]
  group <- group[nna]

  lev  <- levels(group)
  nlev <- length(lev)
  lc   <- p * nlev
  curves <- list()
  clev   <- rep('', lc)  # for each curve what is the level of group

  xmin <- ymin <- 1e30; xmax <- ymax <- -1e30
  ic <- 0
  for(k in 1:p) {
    y <- Y[, k]
    for(g in lev) {
      ic <- ic + 1
      s <- group == g
      z <- switch(method, 
                  lowess=lowess(x[s], y[s], iter=iter, f=f),
                  supsmu=supsmu(x[s], y[s], bass=bass),
                  raw=approx(x[s], y[s], xout=sort(unique(x[s]))))
      
      if(missing(trim))
        trim <- if(sum(s) > 200) 10 / sum(s) else 0
      
      if(trim > 0 && trim < 1) {
        xq <- quantile(x[s], c(trim, 1 - trim))
        s <- z$x >= xq[1] & z$x <= xq[2]
        z <- list(x=z$x[s], y=z$y[s])
      }
      
      if(length(evaluate)) {
        rx   <- range(z$x)
        xseq <- seq(rx[1], rx[2], length=evaluate)
        z <- approx(z, xout=xseq)
      }
      
      if(!missing(fun)) {
        yy <- fun(z$y)
        s <- !is.infinite(yy) & !is.na(yy)
        z <- list(x=z$x[s], y=yy[s])
      }
      
      clev[ic] <- g
      lab <- if(p == 1) g
       else if(nlev == 1 & p == 1) '1'
       else if(nlev == 1 & p > 1) colnames(Y)[k]
       else paste(colnames(Y)[k], g)

      curves[[lab]] <- z
      xmin <- min(xmin, z$x); xmax <- max(xmax, z$x)
      ymin <- min(ymin, z$y); ymax <- max(ymax, z$y)
    }
  }
  if(add) {
    if(missing(xlim))
      xlim <- if(grid) current.panel.limits()$xlim else par('usr')[1:2]
  }
  else {
    if(missing(xlab))
      xlab <- label(x, units=TRUE, plot=TRUE, default=nam[1])
 
    if(missing(xlim)) xlim <- c(xmin, xmax)
    if(missing(ylim)) ylim <- c(ymin, ymax)
    plot(xmin, ymin, xlim=xlim, ylim=ylim,
         type='n', xlab=xlab, ylab=ylab)
  }
  
  lty <- rep(lty, length=lc)
  col <- rep(col, length=lc)
  if(missing(lwd) && is.list(label.curves) && length(label.curves$lwd))
    lwd <- label.curves$lwd
  
  lwd <- rep(lwd, length=lc)

  for(i in 1 : lc) {
    cu <- curves[[i]]
    s <- cu$x >= xlim[1] & cu$x <= xlim[2]
    curves[[i]] <- list(x=cu$x[s], y=cu$y[s])
  }
  if(lines.)
    for(i in 1 : lc)
      gfun$lines(curves[[i]], lty=lty[i], col=col[i], lwd=lwd[i])
  
  if(datadensity) {
    for(i in 1 : nlev) {
      s <- group == lev[i]
      x1 <- x[s]
      for(ii in which(clev == lev[i])) {
        y.x1 <- approx(curves[[ii]], xout=x1)$y
        sopts <- c(list(x=x1, y=y.x1, col=col[ii], grid=grid), scat1d.opts)
        do.call('scat1d', sopts)
      }
    }
  }

  if((is.list(label.curves) || label.curves) && 
     lc > 1 && (!missing(prefix) | !add | !missing(label.curves))) 
    labcurve(curves, lty=lty, col.=col, opts=label.curves, grid=grid)
  
  invisible(curves)
}


panel.plsmo <- function(x, y, subscripts, groups=NULL, type='b', 
                        label.curves=TRUE,
                        lwd  = superpose.line$lwd, 
                        lty  = superpose.line$lty, 
                        pch  = superpose.symbol$pch, 
                        cex  = superpose.symbol$cex, 
                        font = superpose.symbol$font, 
                        col  = NULL, scat1d.opts=NULL, ...)
{
  superpose.symbol <- trellis.par.get("superpose.symbol")
  superpose.line   <- trellis.par.get("superpose.line")
  if(length(groups)) groups <- as.factor(groups)
  
  g  <- unclass(groups)[subscripts]
  ng <- if(length(groups)) max(g) else 1
  
  lty  <- rep(lty, length = ng)
  lwd  <- rep(lwd, length = ng)
  pch  <- rep(pch, length = ng)
  cex  <- rep(cex, length = ng)
  font <- rep(font, length = ng)
  if(!length(col))
    col <- if(type == 'p') superpose.symbol$col else superpose.line$col
  
  col <- rep(col, length = ng)
  lc <-
    if(is.logical(label.curves)) {
      if(label.curves)
        list(lwd=lwd, cex=cex[1])
      else FALSE
    } else c(list(lwd=lwd, cex=cex[1]), label.curves)
  
  if(type != 'p') if(ng > 1)
    plsmo(x, y, group=groups[subscripts, drop=FALSE], 
          add=TRUE, lty=lty, col=col, label.curves=lc, grid=TRUE,
          scat1d.opts=scat1d.opts, ...)
  else
    plsmo(x, y, add=TRUE, lty=lty, col=col, label.curves=lc, grid=TRUE,
          scat1d.opts=scat1d.opts, ...)

  if(type != 'l') {
    if(ng > 1)
      panel.superpose(x, y, subscripts,
                      as.integer(groups),
                      lwd=lwd, lty=lty, pch=pch, cex=cex, 
                      font=font, col=col)
    else
      panel.xyplot(x, y, 
                   lwd=lwd, lty=lty, pch=pch, cex=cex, 
                   font=font, col=col)
    
    if(ng > 1) {
      Key <- function(x=NULL, y=NULL, lev, cex, col, font, pch){
        oldpar <- par(usr=c(0, 1, 0, 1), xpd=NA)
        on.exit(par(oldpar))
        if(is.list(x)) {
          y <- x[[2]]
          x <- x[[1]]
        }
            
        ## Even though par('usr') shows 0,1,0,1 after lattice draws
        ## its plot, it still needs resetting
        if(!length(x))
          x <- 0
        
        if(!length(y))
          y <- 1  ## because of formals()
        
        rlegend(x, y, legend=lev, cex=cex, col=col, pch=pch)
        invisible()
      }
      
      formals(Key) <- list(x=NULL,y=NULL,lev=levels(groups), cex=cex,
                           col=col, font=font, pch=pch)
      .setKey(Key)
    }
  }
}