"prism"<-function(x=matrix(rnorm(120),20,6),...)
{
# prism: creates a (n-1) by (n-1) grid of plots.
# The upper diagonal part comprises scatterplots of columns 2:n of x from left
# to right versus columns 1:(n-1) of x from top to bottom.  The lower diagonal
# part is the PRISM associated with x.
#
# Arguments:
#    x: matrix whose columns are the variables to plot against each other
#  ...: any extra parameters to pass to plotting functions
#       (for example, pch=2)
#
# Written by Jie Li. This function is an integration of "mypairs", "myPairs" written 
# by John Castelloe with some modifications. 



# For testing purposes only (set to F if you don't need "2 vs. 1" labels):
  showlabels <- F

# Make sure x is in matrix form
  x <- as.matrix(x)
  n <- ncol(x)

# Make sure we clean up our mess (revert back to original "par" parameters)
# when we're done
  oldpar <- par("pty", "oma", "mar", "cex", "tck",
    "mgp", "mex", "mfrow")
  on.exit({par(oldpar)})

# Set layout, size and margin parameters
  # make plots square
  par(pty="s")
  # determine new optimal plotting character size
  CEX <- par("cex") * max(7.7/(2*(n-1)+3), 0.6)*0.5
  # specify grid of (n-1) by (n-1) plots
  par(mfcol = c(n-1,n-1))
  # magically cause plots to pack together nicely
  par(oma = rep(3,4))
  par(mar = rep(0.3,4))
  dif <- diff(par("fin"))/2
  if(dif > 0)
    par(omi = c(dif*(n-1), 0, dif*(n-1), 0) + par("omi"))
  else
    par(omi = c(0, (-dif)*(n-1), 0, (-dif)*(n-1)) + par("omi"))
  # re-specify plotting character size (since par(mfrow) screwed it up)
  par(cex = CEX)
  # specify the points type
  par(pch=16)

  
  dat <- x
  m <- n-2
  mat.x <- NULL
  mat.y <- NULL

  for (k in (1:m))
  {
   fres <- sapply(1:(m+1-k), function(i) residuals(lm(dat[,i] ~ dat[,((i+1):(i+k))])))
   bres <- sapply(1:(m+1-k), function(i) residuals(lm(dat[,i+k+1] ~ dat[,((i+1):(i+k))])))
   mat.x<-cbind(mat.x,fres)
   mat.y<-cbind(mat.y,bres)
  }
  
  locmat<-matrix(rep(0,m^2),nrow=m)
  count<-0
  for (i in (1:m))
    for (j in (1:(m+1-i)))
       {
        count<-count+1
        locmat[i,j]<-count
       }

  
# Make plots (i from top to bottom, j from left to right)
  for (i in (2:n))
  {
    for(j in (1:(n-1)))
    {
	  if (j < i)
	  {
	    # prepare an empty plot with correct limits and no axes
		plot(range(x[!is.na(x[,j]),j]), range(x[!is.na(x[,i]),i]),
	       type = "n", axes = F, ...)
		# create empty plot and put border around it
        box()
		# show points for j'th vs. i'th columns of x
	   	points(as.vector(x[,j]),as.vector(x[,i]), ...)
		# show which variables are being plotted
		if (showlabels)
		  text(mean(range(x[!is.na(x[,j]),j])),
		    mean(range(x[!is.na(x[,i]),i])),
		    paste(j,"vs.",i),cex=1.5*CEX)
	  }
	  else
	  {
	  # Get the information of which variables to be plotted against each other on this grid
    loc<-locmat[j-i+1,i-1]
	    # prepare an empty plot with correct limits and no axes
		plot(range(mat.x[!is.na(mat.x[,loc]),loc]), range(mat.y[!is.na(mat.y[,loc]),loc]),
	       type = "n", axes = F, ...)
		# create empty plot and put border around it
    box()
		# show points for j'th vs. i'th columns of x
    points(as.vector(mat.x[,loc]),as.vector(mat.y[,loc]), ...)

	  }
	}
  }

}