# Usage:
# xyaplot (y~x, groups=f, data=df,
#         axes=ay~ax, data.axes=ef, groups.axes=h,
#         axes.spec=list(origin="ao",  min="am", max="aM", interval=TRUE, ticks=TRUE),
#         label="label", label.axes="alabel")
# y~x is the formula for a scatterplot of cases, taken from data, with optional groups
# ay~ax is the formula for the axes, taken from data.axes, with optional groups.axes
# Can specify the origin, min, max of the axis: this means
# (ay,ay) is a unit step in the positive direction;
# min and max say how far to draw the axis for;
# origin says where the planar origin is on this axis
# Note: origin, min, max are characters referring to columns in data.axes
# interval, ticks: either boolean, or referring to columns
# Say whether to draw a line from min to max, and whether to draw tick marks
# label, label.axes label the points and the axes
# style=function(g) returns the style for points in group g
# style.axes=function(g) returns the style for an axis line in group g

equal.num <- function(x,y) identical(all.equal(x, y), TRUE)
reverse <- function(x) x[length(x):1]

biplot.toxy <- function(lambda, origin, dir) dir*(lambda-origin)
biplot.normal <- function(d) { nrm <- c(-d[[2]],d[[1]]); if (nrm[[2]]>=0) nrm else -nrm }
biplot.eucn <- function(d)  sqrt(d[[1]]^2+d[[2]]^2) 
biplot.extent <- function(op,d,xlim,ylim) {
  bd1 <- if (d[[1]]==0) c(Inf,Inf)
    else if (d[[1]]>0) (xlim-op[[1]])/d[[1]]
    else (reverse(xlim)-op[[1]])/d[[1]]
  bd2 <- if (d[[2]]==0) c(Inf,Inf)
    else if (d[[2]]>0) (ylim-op[[2]])/d[[2]]
    else (reverse(ylim)-op[[2]])/d[[2]]
  mM <- c(max(bd1[[1]],bd2[[1]]),min(bd1[[2]],bd2[[2]]))
  mM
  }

princompfull <- function(m) {
  mins <- apply(m,2,min)
  maxs <- apply(m,2,max)
  m <- scale(m, center=TRUE, scale=FALSE)
  means <- attr(m,"scaled:center")
  decomp <- svd(m)
  xhat <- decomp$u[,1:2] %*% diag(decomp$d)[1:2,1:2]
  casedf <- data.frame(label=I(dimnames(m)[[1]]),x=xhat[,1],y=xhat[,2])
  vp <- decomp$v[,1:2]
  u <- apply(vp,1,function(x) sum(x*x))
  vardf <- data.frame(label=I(dimnames(m)[[2]]),x=vp[,1]/u[1],y=vp[,2]/u[2],
                      o=means, min=mins, max=maxs)
  list(observations=casedf, variables=vardf)
}

pointstyle.default <- function(g) {
  n <- length(levels(g))
  i <- which(levels(g)==g)
  sym <- if (n>1) trellis.par.get("superpose.symbol") else trellis.par.get("plot.symbol")
  cex <- rep(sym$cex,length.out=i)[i]
  col <- rep(sym$col,length.out=i)[i]
  font <- rep(sym$font, length.out=i)[i]
  pch <- rep(sym$pch, length.out=i)[i]
  list(cex=cex,col=col,font=font,pch=pch)
}

panel.xyaplot <- function(data, x, y, group, label=NULL, style=pointstyle.default,...) {
  sty <- style(group)
  if (is.null(label))
    panel.xyplot(data[[x]],data[[y]], cex=sty$cex, col=sty$col, pch=sty$pch)
  else
    ltext(data[[x]],data[[y]],data[[label]], cex=sty$cex, col=sty$col, font=sty$font)
}

prepanel.default.xyaplot.axes <- function(data, x,y, axes.spec=list(), ...) {
  n <- nrow(data)
  xn <- data[[x]]
  yn <- data[[y]]
  nmin <- axes.spec$min
  nmax <- axes.spec$max
  norig <- axes.spec$orig
  amin <- if(!is.null(nmin)) data[[nmin]] else NULL
  amax <- if(!is.null(nmax)) data[[nmax]] else NULL
  aorig <- if(!is.null(norig)) data[[norig]] else NULL
  if (is.null(norig) & (!is.null(nmin) | !is.null(nmax))) stop("Must specify origin before bounds")
  xylim <- c(NA,NA,NA,NA)
  ensurepoint <- function(xylim,p) {
    xlim <- range(xylim[1:2],p[[1]], na.rm=TRUE);
    ylim <- range(xylim[3:4],p[[2]], na.rm=TRUE);
    c(xlim,ylim)
    }
  if (!is.null(norig)) for (i in 1:n) {
    d <- c(xn[[i]],yn[[i]])
    o <- aorig[[i]]
    m <- if (!is.null(amin)) amin[[i]] else NA
    M <- if (!is.null(amax)) amax[[i]] else NA
    if (!is.na(o)) {
#      xylim <- ensurepoint(xylim, biplot.toxy(0,o,d))
      if (!is.na(m)) xylim <- ensurepoint(xylim, biplot.toxy(m,o,d))
      if (!is.na(M)) xylim <- ensurepoint(xylim, biplot.toxy(M,o,d))
      }
    }
  res <- list(xlim=xylim[1:2], ylim=xylim[3:4], dx=NA, dy=NA)
  res
}

dupat <- function(v,i) rep(v,length.out=i)[i]

whiten <- function(col,r) {
  c <- col2rgb(col)
  c <- c[,1]
  c <- (1-r)*(255-c)
  c <- 255-c
  rgb(c[1],c[2],c[3],max=255)
}

axesstyle.default <- function(g) {
  n <- length(levels(g))
  i <- which(levels(g)==g)
  sym <- if (n>1) trellis.par.get("superpose.symbol") else trellis.par.get("plot.symbol")
  lin <- if (n>1) trellis.par.get("superpose.line") else trellis.par.get("add.line")
  labelstyle <- list(cex=dupat(sym$cex,i), col=dupat(sym$col,i), font=dupat(sym$font,i))
  intstyle <- list(col=dupat(lin$col,i), lty=dupat(lin$lty,i), lwd=dupat(lin$lwd,i))
  axisstyle <- intstyle
  axisstyle$col <- whiten(axisstyle$col,.8)
  axlabstyle <- trellis.par.get("axis.text")
  axlabstyle$cex <- axlabstyle$cex*0.75
  list(label=labelstyle, axis=axisstyle, interval=intstyle, axislabel=axlabstyle)
}


panel.xyaplot.axes <- function(data, x,y,group, label.axes=NULL, axes.spec=list(), eps=.01,
                               style.axes=axesstyle.default,...) {
  # Work out dimensions
  xlim <- current.viewport()$xscale
  ylim <- current.viewport()$yscale
  xlims <- mean(xlim)+c(-1,1)/2*0.8*(xlim[[2]]-xlim[[1]])
  ylims <- mean(ylim)+c(-1,1)/2*0.8*(ylim[[2]]-ylim[[1]])
  eps <- eps*max(xlim[[2]]-xlim[[1]],ylim[[2]]-ylim[[1]])
  # Basic data and sanity checks
  xn <- data[[x]]
  yn <- data[[y]]
  n <- nrow(data)
  label.axes <- if (!is.null(label.axes)) data[[label.axes]]
  style <- style.axes(group)
  norig <- axes.spec$origin
  nmin <- axes.spec$min
  nmax <- axes.spec$max
  amin <- if (!is.null(nmin)) data[[nmin]] else NULL
  amax <- if (!is.null(nmax)) data[[nmax]] else NULL
  aorig<- if (!is.null(norig)) data[[norig]] else NULL
  interval <- axes.spec$interval
  ticks <- axes.spec$ticks
  if (is.null(interval)) interval <- TRUE
  if (is.null(ticks)) ticks <- TRUE
  if (is.logical(interval)) interval <- rep(interval,length.out=n)
  if (is.logical(ticks)) ticks <- rep(ticks,length.out=n)
  if (is.character(interval)) interval <- data[[interval]]
  if (is.character(ticks)) ticks <- data[[ticks]]
  if (is.null(aorig) & (!is.null(amin) | !is.null(amax))) stop("Must specify origin before bounds")
  for (i in 1:n) {
    # Decide if we're drawing ticks and interval
    tic <- ticks[[i]]
    int <- interval[[i]]
    # Work out the significant points
    d <- c(xn[[i]],yn[[i]])
    epsd <- eps*biplot.normal(d)/biplot.eucn(biplot.normal(d))
    o <- if (is.null(aorig)) NA else aorig[[i]]
    op <- if (is.na(o)) NULL else biplot.toxy(0,o,d)
    m <- if (is.null(amin)) NA else amin[[i]]
    M <- if (is.null(amax)) NA else amax[[i]]
    lab <- if (is.null(label.axes)) NA else label.axes[[i]]
    haveint <- (!is.na(o) & !is.na(m) & !is.na(M))
    axisstyle <- if (int & haveint) style$axis else style$interval
    # Draw the basic axis
    panel.abline(0,d[[2]]/d[[1]],
                 col=axisstyle$col, lty=axisstyle$lty, lwd=axisstyle$lwd)
    # If we want to, and can, draw the interval
    if (int & haveint)
      lline(biplot.toxy(m,o,d),biplot.toxy(M,o,d),
            col=style$interval$col, lty=style$interval$lty, lwd=style$interval$lwd)
    # If we want ticks, draw them
    tickstyle <- style$interval
    tmM <- if (tic & haveint) c(m,M) else biplot.extent(op,d,xlim,ylim)
    tm <- tmM[[1]]
    tM <- tmM[[2]]
    if (tic & !equal.num(tm,tM)) {
      ticksn <- ticks.default(tm,tM)
      adj <- ticktext.direction(d[[1]],d[[2]])
      for (t in ticksn) {
        tp <- biplot.toxy(t,o,d)
        lline(tp-epsd,tp,
              col=tickstyle$col, lty=tickstyle$lty, lwd=tickstyle$lwd)
        dtext(tp[[1]]-epsd[[1]],tp[[2]]-epsd[[2]],as.character(t), adj=adj,
              col=style$axislabel$col, cex=style$axislabel$cex, font=style$axislabel$font)
        }
      }
    # Now draw the labels: at M, or at extent of viewport
    if (!is.na(lab)) {
      labanchor <- NULL
      if (!is.na(o) & !is.na(M)) labanchor <- biplot.toxy(M,o,d)
      if (is.null(labanchor)) {
        op <- biplot.toxy(0,0,d)
        M <- biplot.extent(op,d, xlims, ylims)[[2]]
        labanchor <- biplot.toxy(M,0,d)
        }
      th <- atan2(d[[2]],d[[1]])*360/(2*pi)
      if (d[[1]]<0) th <- th-180
      dtext(labanchor[[1]]+epsd[[1]],labanchor[[2]]+epsd[[2]],lab,adj=c(.5,0),srt=th,
            cex=style$label$cex, col=style$label$col, font=style$label$font)
      }
    }
  }

rotatepoint <- function(p,th)
  (p %*% matrix(c(cos(th),-sin(th),sin(th),cos(th)),ncol=2))[1,]


dtext <- function (x, y = NULL, labels = seq(along = x), col = add.text$col, 
    cex = add.text$cex, srt = 0, font = 1, adj = 0.5, ...) {
    add.text <- trellis.par.get("add.text")
    adj <- rep(adj, length.out=2)
    adjx <- adj[[1]]
    adjy <- adj[[2]]
    justx <- if (adjx==0) "left" else if (adjx==1) "right" else "centre"
    justy <- if (adjy==0) "bottom" else if (adjy==1) "top" else "centre"
    just <- c(justx,justy)
    xy <- xy.coords(x, y)
    grid.text(label = labels, x = xy$x, y = xy$y, gp = gpar(col = col, 
        font = font, fontsize = cex * trellis.par.get("fontsize")$default), 
        just = just, rot = srt, default.units = "native")
}

ticktext.direction <- function(dx,dy) {
  if (dy<0) { dx <- -dx; dy <- -dy }
  th <- atan2(dy,dx)
  if (equal.num(th,pi/2))
    c(1,.5)
  else if (equal.num(th,0) || equal.num(th,pi))
    c(.5,1)
  else if (th<pi/2)
    c(0,1)
  else
    c(1,1)
}


lline <- function(pt1,pt2,...)
  llines(c(pt1[[1]],pt2[[1]]),c(pt1[[2]],pt2[[2]]),...)

ticks.default <- function(m,M) {
  d <- sign(M+m)
  r <- M-m;
  i <- 10^floor(log10(r))
  n <- ticks.visible(m,M,i)
  if (n[[2]]-n[[1]]>5) i <- 2*i
  if (n[[2]]-n[[1]]<4) i <- i/2
  if (n[[2]]-n[[1]]<2) i <- 2*i/5
  n <- ticks.visible(m,M,i)
  if (n[[2]]-n[[1]]>5) i <- 2*i
  n <- ticks.visible(m,M,i)
  n[[1]]:n[[2]] * i
}

ticks.visible <- function(m,M,i) {
  n0 <- if(m<0) -tick.before(-m,i) else tick.after(m,i)
  n1 <- if(M<0) -tick.after(-M,i) else tick.before(M,i)
  c(n0,n1)
}

tick.before <- function(x,i) floor(x/i)
tick.after <- function(x,i) ceiling(x/i)

panel.superpose.ex <- function(data, groups, panel.groups, ...) {
  gru <- sort(unique(groups))
  ngr <- length(gru)
  for (i in seq(along=gru)) {
    gr <- gru[i]
    sel <- (groups==gr)
    if (any(sel)) {
      args <- list(data=data[sel,], group=gr, ...)
      do.call("panel.groups", args)
      }
    }
  }

panel2.default <- function(data, groups, panel, x, y,
                   data.axes, groups.axes, panel.axes, x.axes, y.axes,...) {
  panel.axes(data=data.axes,groups=groups.axes,x=x.axes,y=y.axes,...)
  panel(data=data,groups=groups,x=x,y=y,...)
  }

prepanel2.default <- function(data, x,y, prepanel,
                      data.axes, x.axes, y.axes, prepanel.axes, ...) {
  prep1 <- prepanel(data[[x]],data[[y]],...)
  prep2 <- prepanel.axes(data=data.axes,x=x.axes,y=y.axes,...)
  xlim <- range(prep1$xlim, prep2$xlim, na.rm=TRUE)
  ylim <- range(prep1$ylim, prep2$ylim, na.rm=TRUE)
  p <- list(xlim=xlim, ylim=ylim, dx=prep1$dx, dy=prep1$dy)
  p
  }

as.reallyFactorOrShingle <- function(x,...) {
  y <- as.factorOrShingle(x,...)
  if (!is.factor(y) & !is.shingle(y)) y <- as.factor(y)
  y
}

xyaplot <- function (formula, axes, groups=1, groups.axes=1, data, data.axes,
                     subset=TRUE, subset.axes=TRUE,
                     panel2=panel2.default, prepanel2=prepanel2.default,
                     panel=function(...) panel.superpose.ex(...,panel.groups=panel.xyaplot),
                     panel.axes=function(...) panel.superpose.ex(...,panel.groups=panel.xyaplot.axes),
                     prepanel=prepanel.default.xyplot, prepanel.axes=prepanel.default.xyaplot.axes,
                     aspectfunction=default.aspectfunction,
                     aspect = "fill", layout = NULL,
                     scales=list(draw=FALSE), strip=function(...) strip.default(...,strip.names=TRUE),
                     xlab="",xlim,ylab="",ylim,...) {
  dots <- list(...)
  # Check sanity of groups and subset arguments
  ndata <- nrow(data)
  naxes <- nrow(data.axes)
  groups <- eval(substitute(groups), data, parent.frame())
  groups.axes <- eval(substitute(groups.axes), data.axes, parent.frame())
  groups <- as.factor(groups)
  groups.axes <- as.factor(groups.axes)
  if (length(groups)==1) groups <- rep(groups, length.out=ndata)
  if (length(groups.axes)==1) groups.axes <- rep(groups.axes, length.out=naxes)
  if (length(subset)==1) subset <- rep(subset, length.out=ndata)
  if (length(subset.axes)==1) subset.axes <- rep(subset.axes, length.out=naxes)
  if (length(groups)!=ndata) stop("Groups argument of wrong length")
  if (length(groups.axes)!=naxes) stop("Groups.axes argument of wrong length")
  if (length(subset)!=ndata) stop("Subset argument of wrong length")
  if (length(subset.axes)!=naxes) stop("Subset.axes argument of wrong length")
  # Something about the labels being lists?
  # Look at data and formulae and conditioning
  form <- latticeParseFormula(formula, data)
  form.axes <- latticeParseFormula(axes, data.axes)
  x <- form$left.name
  y <- form$right.name
  x.axes <- form.axes$left.name
  y.axes <- form.axes$right.name
  # Work out the layout of the panels
  cond <- lapply(form$condition, as.reallyFactorOrShingle, subset, drop=TRUE)
  cond.axes <- lapply(form.axes$condition, as.reallyFactorOrShingle, subset.axes, drop=TRUE)
  condn <- length(cond)
  condn.axes <- length(cond.axes)
  cond.all <- c(cond,cond.axes)
  condn.all <- length(cond.all)
  if (condn.all==0) strip <- FALSE
  cond.maxlevel <- if (condn.all>0) unlist(lapply(cond.all, nlevels)) else 1
  layout <- compute.layout(layout, cond.maxlevel, skip=FALSE)
  plots.per.page <- max(layout[1] * layout[2], layout[2]) # since perhaps layout[1]=0
  number.of.pages <- layout[3]
  nplots <- plots.per.page * number.of.pages
  # Marshal common arguments
  tskargs <- c(list(aspect=aspect, strip=strip, panel=panel2, xlab=xlab, ylab=ylab), dots)
  foo <- do.call("trellis.skeleton",tskargs)
  dots <- foo$dots
  foo <- foo$foo
  foo <- c(foo, do.call("construct.scales", scales))
  foo$call <- match.call()
  foo$condlevels <- if (condn.all==0) list("alllevels") else lapply(cond.all, levels)
  foo$skip <- rep(FALSE, length=plots.per.page)
  foo$layout <-  layout
  panel2args <- list(panel=panel, panel.axes=panel.axes)
  foo$panel.args.common <- c(list(x=x,y=y,x.axes=x.axes,y.axes=y.axes),panel2args,dots)
  # Marshal arguments for each of the different panels
  foo$panel.args <- as.list(1:nplots)
  cond.level <- rep(1, max(1,condn.all))
  panel.number <- 1
  for (page.number in 1:number.of.pages)
    if (all(cond.level<=cond.maxlevel))
      for (plot in 1:plots.per.page) {
        if (foo$skip[plot]) foo$panel.args[[panel.number]] <- FALSE
        else if (all(cond.level<=cond.maxlevel)) {
          # Select the indices which will be displayed in this panel (specified by cond.level)
          id <- subset
          if (condn>0) for (i in 1:condn) {
            var <- cond[[i]]
            idr <- if (is.shingle(var))
              (var>=levels(var)[[cond.level[i]]][1] & var<=levels(var)[[cond.level[i]]][2])
            else (as.numeric(var)==cond.level[i])
            id <- id & idr
            }
           id.axes <- subset.axes
          if (condn.axes>0) for (i in 1:condn.axes) {
            var <- cond.axes[[i]]
            idr <- if (is.shingle(var))
              (var>=levels(var)[[cond.level[condn+i]]][1] & var<=levels(var)[[cond.level[condn+i]]][2])
            else (as.numeric(var)==cond.level[condn+i])
            id.axes <- id.axes & idr
            }
          panel.data <- data[id,]
          panel.data.axes <- data.axes[id.axes,]
          panel.groups <- groups[id]
          panel.groups.axes <- groups.axes[id.axes]
          panel.args <- list(data=panel.data, data.axes=panel.data.axes,
                             groups=panel.groups, groups.axes=panel.groups.axes)
          foo$panel.args[[panel.number]] <- panel.args
          cond.level <- cupdate(cond.level, cond.maxlevel)
          }
        panel.number <- panel.number+1
      }
  prep2 <- function(...) prepanel2(...,prepanel=prepanel, prepanel.axes=prepanel.axes)
  lims <- limits.and.aspect(prep2, prepanel=NULL,
                            xlim=xlim, ylim=ylim, have.xlim=!missing(xlim), have.ylim=!missing(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)
  res <- aspectfunction(aspect,lims$x.limits,lims$y.limits)
  foo$aspect.fill <- FALSE
  lims$aspect.ratio <- res$aspect
  lims$x.limits <- res$xlim
  lims$y.limits <- res$ylim
  foo <- c(foo,lims)
  class(foo) <- "trellis"
  foo
}

default.aspectfunction <- function(aspect,xlim,ylim) {
  dx <- xlim[[2]]-xlim[[1]]
  dy <- ylim[[2]]-ylim[[1]]
  if (aspect=="fill")
    list(aspect=dy/dx,xlim=xlim,ylim=ylim)
  else {
    d <- max(dx,dy/aspect)
    xlim <- mean(xlim)+c(-d,d)/2
    ylim <- mean(ylim)+aspect*c(-d,d)/2
    list(aspect=aspect,xlim=xlim,ylim=ylim)
  }
}