#$Author: sinnwell $
#$Date: 2010/01/08 19:48:53 $
#$Header: /people/biostat3/sinnwell/Projects/LDpairs/Make/RCS/plot.ld.pairs.q,v 1.10 2010/01/08 19:48:53 sinnwell Exp $
#$Locker:  $
#$Log: plot.ld.pairs.q,v $
#Revision 1.10  2010/01/08 19:48:53  sinnwell
#put a title of which measure used, check if r2 is null before using it
#
#Revision 1.9  2009/12/28 15:26:48  sinnwell
#fix oldpar to work for R and S+
#
#Revision 1.8  2009/12/23 22:29:34  sinnwell
#version to work in both R and Splus
#
#Revision 1.7  2007/02/23 14:38:52  sinnwell
#plot either dprime.ave, r2, or -log10(comp.pval)
#only 5 categories for each
#
#Revision 1.6  2004/09/01 19:11:36  sinnwell
#change default plot size to figure min.pct
#
#Revision 1.5  2004/05/25 20:29:35  sinnwell
#add srt and adj for options with locus.labels
#
#Revision 1.4  2004/05/21 20:21:37  sinnwell
#fix parameter checks, make option for locus.label on axes
#
#Revision 1.3  2004/03/30 17:00:55  sinnwell
#change upper range from 1.0 to 1.1
#
#Revision 1.2  2004/03/25 22:35:28  sinnwell
#add parameters fill.type and line.density
#
#Revision 1.1  2003/10/14 19:51:04  sinnwell
#Initial revision
#
plot.ld.pairs <- function(x, ld.measure="dprime.ave", line.density=6, fill.type="both",
        locus.label=as.character(1:max(c(x$loc1,x$loc2))), adj=.5, srt=0, ...)

# Schaid, DJ and Sinnwell, JP
# Mayo Division of Biostatistics
# 5/2003

# updated to plot dprime.ave, r2, or -log10(comp.pval)
# with 5 categories
# JP Sinnwell 2/2007
{
  # control density of lines with line.density
  # contrast*line.density = number of lines difference in a 10% ld level
  # change per 1 square inch of the box

  # dots <- as.list(substitute(list(...)))[-1]
  # nm <-  names(dots)
  oldpar <- par()
  if(is.R()) {
    oldpar$cin <- oldpar$cra <- oldpar$csi <- oldpar$cxy <- oldpar$din <- oldpar$gamma <- NULL
  }
  on.exit(par(oldpar))
  ld.meas.options = c("dprime.ave", "r2", "comp.pval")
  ld.meas.num <- pmatch(ld.measure, c("dprime.ave", "r2", "comp.pval"))
  if(is.na(ld.meas.num)) {
    warning(paste("Unrecognized ld.measure: ", ld.measure, ", dprime.ave is used.\n"))
    ld.meas.num <- 1
  }
  if(ld.meas.num==2 & is.null(x$r2)) {
    warning("r2 not calculated, so default of dprime.ave is plotted\n")
    ld.meas.num <- 1
  }

  ld.meas.used = ld.meas.options[ld.meas.num]
  
  if(line.density < 1) {
    cat("Warning: line.density out of range, using default 6 \n")
    line.density=6
  }

  fill.type <- pmatch(fill.type,c("up","down","both"))
  if(is.na(fill.type)) {
    cat("Warning: fill.type should be up, down, or both; using both \n")
    fill.type=3
  }

  # function to plot the outline of a box, given coordinates of the center
  plot.box <- function(x.center, y.center, width, height=width){
     w <- width/2
     h <- height/2
     x1 <- x.center - w
     x2 <- x.center + w
     y1 <- y.center - h
     y2 <- y.center + h
     lines(c(x1,x2),c(y1,y1))
     lines(c(x1,x2),c(y2,y2))
     lines(c(x1,x1),c(y1,y2))
     lines(c(x2,x2),c(y1,y2))
   }

  # function to fill box from plot.box with diagonal lines.
  fill.box <- function(y.center, x.center, width, density, fill.type){

    # density is number of lines to be put in each half divided by the diagonals
    # default is to plot on a square block, but with h!=w, it can plot rectangles
    # type is lines directed either up, down, or both at the same time 

     w <- width/2
     x1 <- x.center - w
     x2 <- x.center + w
     y1 <- y.center - w
     y2 <- y.center + w
     delta <- 1/density

     d <- ceiling(density)
     for(i in 1:d) {

       if(fill.type==1 | fill.type==3) {
         # upper left triangle, slope up
         y1.new <- y1 + (i-1)*delta
         x2.new <- x1 + y2 - y1.new
         if(x2.new > x1)
           lines(c(x1,x2.new),c(y1.new,y2))

         # lower right triangle, slope up
         if(i>=1){
           y2.new <- y2 - (i-1)*delta
           x1.new <- x2 - y2.new + y1
           if(x1.new <  x2)
             lines(c(x1.new,x2),c(y1,y2.new))
         }
       }

       if(fill.type==2 | fill.type==3) {
          # upper right triange, slope down
         y1.new <- y1 + (i-1)*delta
         x1.new <- x2 - y2 + y1.new
         if(y1.new< y2)
           lines(c(x1.new,x2),c(y2,y1.new))

           # lower left triangle, slope down
         if(i >= 1){
           y2.new <-  y2 - (i-1)*delta
           x2.new <- x1 + y2.new - y1
           if(y2.new > y1)
             lines(c(x1,x2.new),c(y2.new,y1))
         }
       }
     }
   }

  # find the smallest percentage of a 'standard' plot the current plot is
  # here, call 'standard' plot size pin = 7.132" x 5.032" !!now 7.256 x 5.216!! jps--9/04

   if(is.R()) {
    fullpin <- c(5.754751,5.15524) 
  } else {
    fullpin <- c(7.256, 5.216)
  }
  
  min.pct <- min(oldpar$pin[1]/fullpin[1], oldpar$pin[2]/fullpin[2])
  

  ### PLOT LD PAIRS MATRIX ###
  
  n.loci <- max(c(x$loc1,x$loc2))

  if(length(locus.label) != n.loci) {
    cat("Warning: locus.label not correct length! using", 1:n.loci, "\n")
    locus.label <- 1:n.loci
  }
  y.range <-  c(0.5,n.loci)

  pinLeft <- c(4,4)
  maiLeft <- c(1,1,1,2)
    
  # since plotting in inches, consider oldpar as plotting area to work within
  par(pin=pinLeft*min.pct,mai=maiLeft*min.pct, new=FALSE)
  plot(c(0.5,n.loci),y.range,type="n",xlab="locus",ylab="locus", xaxt="n",yaxt="n")
  axis(1,1:n.loci,locus.label,cex=min.pct, srt=srt, adj=adj)
  axis(2,n.loci:1,locus.label,cex=min.pct, adj=adj)  # tck=0
  title(ld.meas.used)
  
  # height of box for plot, shrink% * 4 inches / no.-of-boxes
  ht1 <- min.pct*4/(n.loci-.5)
  
  # Will plot number of lines per inch of the plot.  Calculate inches in the box.
  # Find half of the diagonal of box in inches
  # use this to decide how many lines to have per half of box
  diag1 <- sqrt(2*(ht1^2))/2
  base.lines <- line.density*diag1

  # categorize ld measure
  # dprime.ave from 0 to .2, .4, etc.  get 0 by starting negative
  # r2 from 0, .2, .4, .6, .8, 1
  # -log10(comp.pval) from 0,.5, 1, 2, 4, 4+
  ld.meas.cat <- switch(ld.meas.num,
                        as.numeric(cut(abs(x$dprime.ave), breaks=c(-1.1,.2,.4,.6,.8,1.1))),
                        as.numeric(cut(x$r2, breaks=c(-.1,.2,.4,.6,.8,1.1))),
                        as.numeric(cut(-log10(x$comp.pval+1e-8), breaks=c(-.01,.5,1,2,4,9))))
                                                   # 1e-8 added for pval of zero
  for(k in 1:length(x$loc1)){
    i <- x$loc1[k]
    j <- x$loc2[k]
    n.lines <- ld.meas.cat[k]*base.lines
    plot.box(y.center=(n.loci - j + 1), x.center=i, width=1)
    fill.box(y.center=(n.loci - j + 1), x.center=i, width=1,
             density = n.lines, fill.type=fill.type)
  }

  ##########  make LEGEND  ################
  label <- switch(ld.meas.num,
                  c("0.0-0.2", "0.2-0.4", "0.4-0.6", "0.6-0.8", "0.8-1.0"),
                  c("0.0-0.2", "0.2-0.4", "0.4-0.6", "0.6-0.8", "0.8-1.0"),
                  c("0-.5", ".5-1", "1-2", "2-4", "4+"))

  y <-  seq(.5, 4.5, by=1)

  # plot legend on the right, measure in inches
  # min.pct is percent to shrink if not a full size plot
  ## put w/ above plot by new=T
  
  pinRight <- c(.5,4)
  maiRight <- c(1,5.5,1,1)
  par(pin=pinRight*min.pct, mai=maiRight*min.pct, new=T, cex=min.pct*.8, new=TRUE)
##  par(pin=c(.4,4)*min.pct,mai=c(1,4,1,0)*min.pct,
##      new=T, cex=min.pct*.8)
  plot(c(0,1),c(0,10), type="n",bty="n", xlab="",ylab="",xaxt="n",yaxt="n")
  axis(2,at=y, labels=label,pos=(-.1)*min.pct)


  # height of boxes, default ht is .4, reduce by min.pct
  ht2 <- min.pct*0.4
  # length of half the diagonal
  diag2 <- sqrt(2*(ht2^2))/2
  base.lines <- diag2*line.density
  for(j in 1:5) {
    n.lines <- j*base.lines
    plot.box(y.center=y[j], x.center=0.5, width=1)
    fill.box(y.center=y[j], x.center=0.5, width=1,
             density=n.lines, fill.type=fill.type)
  }
  
  invisible(x)

}
