#----------------------------------------------------------------------------
#                     BLiP function
#                     version 1.0    July 20, 1994
# updated 9/27/94: add freqtext so that when freq=0, it won't write "0"
#---------------------------------------------------------------------------
blip_function(...,std.plot=NA,graph.region=c(.25,.75),graph.type="c",
graph.width="f",graph.uniform=T,nclass=5,
box.p=c(.25,.5,.75),box.p.bar=rep(T,length(box.p)),
box.label.on=T,box.label.cex=par("cex"),box.label.above=.05,
box.density=0,box.angle=45, box.col=par("col"),
line.p=NA,line.sd=NA,line.se=NA,mean.pch=NA,mean.cex=par("cex"),
mean.col=par("col"),line.col=par("col"),line.width=par("lwd"),
point.pattern="on-line",point.type="p",point.cex=.7,point.pch=par("pch"),
point.col=par("col"),xlim=NA,axes=T,xaxt="s",xlab="",ylab=NA,breaks=NA,
seed=500)
{
#----------------------------------------------------------------------------
#  Extract names of vectors for labeling
#  at y-axis. The by-product is to strip away the extra layer of "list" from
#  the input vector when user enter list(x1,x2) or split(x1,x2).
#  It returns "ylabel"(which is a vector) and "x"(which is a list)
#----------------------------------------------------------------------------
 x_list(...)
 if (is.list(x[[1]])==F)                      #user entered blip(x1,x2)
 {
  if (all(is.na(ylab))==T)
  {
   if (is.null(names(match.call())))         #didn't enter optional arguments
   {
    ylab_as.character(match.call())
    ylab_ylab[ylab!="blip"]
   }
   else                                     #entered optional arguments
   {
    temp_match.call()
    ylab_paste(temp[names(temp)==""])
    ylab_ylab[ylab!="blip"]
   }
  }
 }
 else
 {
  if (is.null(names(x[[1]]))==T)               #entered blip(list(x1,x2))
  {
   if (all(is.na(ylab))==T)
   {
    ylab_vector()
    for (i in 1:length(x[[1]]))
    {
     ylab[i]_i
    }
   }
   temp_list()
   for (i in 1:length(x[[1]]))
   {
    temp[[i]]_x[[1]][[i]]
   }
   x_temp
  }
  else                                        #entered 'split(x1,x2)'
  {
   x_x[[1]]
   if (all(is.na(ylab))==T) ylab_names(x)
  }
 }
 if (length(ylab)!=length(x)) stop("Number of elements in ylab is incorrect")
#---------------------------------------------------------------------------
  minmax_check.common.input(x,std.plot,graph.region,nclass,xlim,breaks)
#----------------------------Define some common variables-------------------
  def_define.var(x,graph.region)
  x_def$x
  nvectors_def$nvectors
  n_def$n
  boxwidth_def$boxwidth
  boxcenter_def$boxcenter
  lower.edge_def$lower.edge
  upper.edge_def$upper.edge
  low_def$low
  up_def$up
  x.min_minmax$x.min
  x.max_minmax$x.max
#--------------------------------------------------------------------------
 if (std.plot=="boxplot")
 {
  bp(x,nvectors,lower.edge,upper.edge,xlim,line.col,line.width,
     box.density,box.angle,box.col,point.cex,point.pch,point.col,
     axes,xaxt,xlab)
 }
 else if (std.plot=="histogram")
 {
  histo(x,nvectors,graph.uniform,low,boxwidth,nclass,
        box.label.cex,box.label.above,box.label.on,xlim,line.col,line.width,
         box.density,box.angle,box.col,axes,xaxt,xlab,breaks,x.min,x.max)
 }
 else if (std.plot=="polygon")
 {
  polygn(x,nvectors,graph.uniform,low,boxwidth,nclass,
         box.label.cex,box.label.above,box.label.on,xlim,line.col,line.width,
         box.density,box.angle,box.col,axes,xaxt,xlab,breaks,x.min,x.max)
 }
 else
 {
  ans1_check.input.blip(x,nvectors,box.p,box.p.bar,graph.type,graph.width,
                        point.pattern,line.p,line.sd,line.se)
  xpoints_ans1$xpoints
  box.bounds_ans1$box.bounds
  xpoints_frame.x(x,nvectors,n,mean.pch,mean.cex,mean.col,
                  line.sd,line.se,xlim,graph.type,
                  boxcenter,low,xpoints,line.width,line.col,
                  axes,xaxt,xlab)
  ans2_ draw.box(x,box.p,box.p.bar,box.bounds,nvectors,n,
                 graph.width,graph.type,graph.uniform,
                 point.pattern,nclass,boxcenter,boxwidth,
                 low,up,box.label.on,box.label.cex,box.label.above,
                 box.density,box.angle,box.col,line.width,line.col)
  draw.line(x,nvectors,line.p,graph.type,boxcenter,low,line.width,line.col)
  draw.point(xpoints,nvectors,n,graph.width,graph.type,graph.uniform,
             point.pattern,ans2,boxwidth,boxcenter,low,up,
             seed,line.width,line.col,point.type,point.cex,point.pch,point.col)
 }
 axis(2,at=boxcenter,labels=ylab,adj=1,ticks=F)
}
 

#-------------------------Include all the functions------------------------
 

check.common.input_function(x,std.plot,graph.region,nclass,xlim,breaks)
{
#--------------------------------------------------------------------------
#           Check input error and put out error message
#--------------------------------------------------------------------------
 if (std.plot != "boxplot" &  std.plot !="histogram" & std.plot !="polygon" &
    !is.na(std.plot))
  stop("Incorrect value for std.plot option")
 if (graph.region[1]>=graph.region[2])
  stop("Value of the lower edge of a box has to be smaller than that of
        upper edge")
 if (length(graph.region) !=2)
  stop("Incorrect values for graph.region option")
 if (is.numeric(nclass)==F)
  stop("Non-numeric value for nclass")
 if (any(!is.na(xlim))==T)
 {
  if (xlim[2]<=xlim[1] | length(xlim)!=2)
   stop("Incorrect value for xlim")
 }
 x.min_min(sapply(x,min))
 x.max_max(sapply(x,max))
 if (std.plot=="histogram" | std.plot=="polygon")
 {
  if (all(!is.na(breaks)))
  {
   if (breaks[1]>x.min | breaks[length(breaks)]<x.max)
    stop("Break points should cover the range of the data")
   if (all(breaks == sort(breaks))==F)
    stop("Break points are not in ascending order")
   x.min_breaks[1]
   x.max_breaks[length(breaks)]
  }
 }
 return(x.min=x.min,x.max=x.max)
}

#---------------------------------------------------------------------------
define.var_function(x,graph.region)
{
#----------------------------Define some common variables------------------
 x_lapply(x,sort)
 nvectors_length(x)
 n_vector()
 low_vector()
 up_vector()
 boxwidth_vector()
 boxcenter_vector()
 lower.edge_graph.region[1]
 upper.edge_graph.region[2]
 for (k in 1:nvectors)
 {
  n[k]_length(x[[k]])
  low[k]_(k-1)/nvectors+(1/nvectors)*lower.edge
  up[k]_(k-1)/nvectors+(1/nvectors)*upper.edge
  boxwidth[k]_up[k]-low[k]
  boxcenter[k]_low[k]+boxwidth[k]/2
 }
 return(x=x,nvectors=nvectors,n=n,low=low,up=up,boxwidth=boxwidth,
        boxcenter=boxcenter,lower.edge=lower.edge,upper.edge=upper.edge)
}

#---------------------------------------------------------------------------
bp_function(x,nvectors,lower.edge,upper.edge,xlim,line.col,
            line.width,box.density,box.angle,box.col,point.cex,point.pch,
            point.col,axes,xaxt,xlab)
{
#------------------------Reset some variables------------------------------
  box.p_c(.25,.5,.75)
  graph.width_"f"
  graph.type_"c"
  point.pattern_"on-line"
  n_vector()
#------------------------Draw frame-----------------------------------------
 if (all(is.na(xlim))==T)
 {
  x.min_min(sapply(x,min))
  x.max_max(sapply(x,max))
  if (x.min>=0) fake.x.min_.95*x.min
  else fake.x.min_1.05*x.min
  if (x.max>=0) fake.x.max_1.05*x.max
  else fake.x.max_.95*x.max
  fake.x_c(fake.x.min,fake.x.max)
 }
 else fake.x_xlim
 plot(fake.x,c(0,1),type="n",axes=F,ylab="",xlab=xlab)
 if (axes==T)
 {
  axis(1,xaxt=xaxt)
  box()
 }
#---------------------Determine the height of the box-----------------------
 for (k in 1:nvectors)
 {
  n[k]_length(x[[k]])
  low_(k-1)/nvectors+(1/nvectors)*lower.edge
  up_(k-1)/nvectors+(1/nvectors)*upper.edge
  boxwidth_up-low
  boxcenter_low+boxwidth/2
  q_quantile(x[[k]],probs=box.p)
  nq_length(q)
  loweredge_rep(low,nq)
  upperedge_rep(up,nq)
#--------------------Draw the box------------------------------------------
  polygon(c(q,rev(q)),c(upperedge,rev(loweredge)),
          density=box.density,col=box.col,angle=box.angle,lwd=line.width,
          border=F)
  segments(q,loweredge,q,upperedge,lwd=line.width,col=line.col)
  lines(q,loweredge,lwd=line.width,col=line.col)
  lines(q,upperedge,lwd=line.width,col=line.col)
#--------------------------Draw lines-------------------------------------
  qr_quantile(x[[k]],probs=c(.25,.75))
  xlines_c(max((qr[1]-1.5*(qr[2]-qr[1])),min(x[[k]])),qr[1],NA,qr[2],
           min((qr[2]+1.5*(qr[2]-qr[1])),max(x[[k]])))
  nxlines_length(xlines)
  lines(xlines,rep(boxcenter,nxlines),lwd=line.width,col=line.col)
  xpoints_c(x[[k]][x[[k]]<xlines[1]],x[[k]][x[[k]]>xlines[nxlines]])
#--------------------------Draw points-----------------------------------
  points(xpoints,rep(boxcenter,length(xpoints)),pch=point.pch,cex=point.cex,
         col=point.col)
 }
}

#-----------------------------------------------------------------------
histo_function(x,nvectors,graph.uniform,low,boxwidth,nclass,
box.label.cex,box.label.above,box.label.on,xlim,line.col,line.width,
box.density,box.angle,box.col,axes,xaxt,xlab,breaks,x.min,x.max)
{
#------------------------Draw frame-----------------------------------------
 if (all(is.na(xlim))==T)
 {
  if (x.min>=0) fake.x.min_.95*x.min
  else fake.x.min_1.05*x.min
  if (x.max>=0) fake.x.max_1.05*x.max
  else fake.x.max_.95*x.max
  fake.x_c(fake.x.min,fake.x.max)
 }
 else fake.x_xlim
 plot(fake.x,c(0,1),type="n",axes=F,ylab="",xlab=xlab)
 if (axes==T)
 {
  axis(1,xaxt=xaxt)
  box()
 }
#---------------------Determine frequency in each interval--------------
 freq_list()
 freqtext_list()
 line_list()
 qlines_list()
 midpoints_list()
 if (graph.uniform==T)
 {
  if (all(is.na(breaks))) intervals_seq(x.min,x.max,by=(x.max-x.min)/nclass)
  else intervals_breaks
  nintervals_length(intervals)
  for (k in 1:nvectors)
  {
   freq[[k]]_as.vector(table(cut(x[[k]],
                   breaks=c(intervals[1]-.1,intervals[2:(nintervals-1)],
                            intervals[nintervals]+.1))))
   freqtext[[k]] _ freq[[k]]
   freqtext[[k]][freq[[k]]==0] _ NA
  }
  max.freq_max(sapply(freq,max))
  for (k in 1:nvectors)
  {
   line[[k]]_matrix(0,(nintervals-1),2)
   line[[k]][,1]_intervals[1:(nintervals-1)]
   line[[k]][,2]_intervals[2:nintervals]
   midpoints[[k]]_line[[k]][,1]+(line[[k]][,2]-line[[k]][,1])/2
   upperedge_low[k]+(freq[[k]]/max.freq)*boxwidth[k]
   loweredge_rep(low[k],2)
   for (i in 1:(nintervals-1))
   {
    polygon(c(line[[k]][i,1],line[[k]][i,1],line[[k]][i,2],line[[k]][i,2]),
            c(loweredge[1],upperedge[i],upperedge[i],loweredge[1]),border=F,
            density=box.density,angle=box.angle,col=box.col,lwd=line.width)
    segments(line[[k]][i,1],loweredge[1],line[[k]][i,1],upperedge[i],
             lwd=line.width,col=line.col)
    segments(line[[k]][i,2],loweredge[1],line[[k]][i,2],upperedge[i],
             lwd=line.width,col=line.col)
    lines(line[[k]][i,1:2],c(upperedge[i],upperedge[i]),
          lwd=line.width,col=line.col)
   }
   lines(c(x.min,x.max),loweredge,lwd=line.width,col=line.col)
   if (box.label.on==T)
   text(midpoints[[k]],upperedge+box.label.above,freqtext[[k]],cex=box.label.cex)
  }
 }
 else
 {
  for (k in 1:nvectors)
  {
   if (all(is.na(breaks))) intervals_seq(min(x[[k]]),max(x[[k]]),
                                     by=(max(x[[k]])-min(x[[k]]))/nclass)
   else intervals_breaks
   nintervals_length(intervals)
   freq[[k]]_as.vector(table(cut(x[[k]],
               breaks=c(intervals[1]-.1,intervals[2:(nintervals-1)],
                        intervals[nintervals]+.1))))
   line[[k]]_matrix(0,(nintervals-1),2)
   line[[k]][,1]_intervals[1:(nintervals-1)]
   line[[k]][,2]_intervals[2:nintervals]
   midpoints[[k]]_line[[k]][,1]+(line[[k]][,2]-line[[k]][,1])/2
   upperedge_low[k]+(freq[[k]]/max(freq[[k]]))*boxwidth[k]
   loweredge_rep(low[k],2)
   for (i in 1:(nintervals-1))
   {
    polygon(c(line[[k]][i,1],line[[k]][i,1],line[[k]][i,2],line[[k]][i,2]),
            c(loweredge[1],upperedge[i],upperedge[i],loweredge[1]),border=F,
            density=box.density,angle=box.angle,col=box.col,lwd=line.width)
    segments(line[[k]][i,1],loweredge[1],line[[k]][i,1],upperedge[i],
             lwd=line.width,col=line.col)
    segments(line[[k]][i,2],loweredge[1],line[[k]][i,2],upperedge[i],
             lwd=line.width,col=line.col)
    lines(line[[k]][i,1:2],c(upperedge[i],upperedge[i]),
          lwd=line.width,col=line.col)
   }
   lines(c(min(x[[k]]),max(x[[k]])),c(loweredge[1],loweredge[1]),
         lwd=line.width,col=line.col)
   if (box.label.on==T)
    text(midpoints[[k]],upperedge+box.label.above,freqtext[[k]],cex=box.label.cex)
  }
 }
}

#---------------------------------------------------------------------------
polygn_function(x,nvectors,graph.uniform,low,boxwidth,nclass,
box.label.cex,box.label.above,box.label.on,xlim,line.col,line.width,
box.density,box.angle,box.col,axes,xaxt,xlab,breaks,x.min,x.max)
{
#------Determine frequency in each interval and draw polygon-----------------
  freq_list()
  line_list()
  qlines_list()
  midpoints_list()
  xpoints_list()
  ypoints_list()
  loweredge_list()
  upperedge_list()
  if (graph.uniform==T)
  {
   x.min_min(sapply(x,min))
   x.max_max(sapply(x,max))
   if (all(is.na(breaks))) intervals_seq(x.min,x.max,by=(x.max-x.min)/nclass)
   else intervals_breaks
   nintervals_length(intervals)
   for (k in 1:nvectors)
   {
    freq[[k]]_as.vector(table(cut(x[[k]],
                breaks=c(intervals[1]-.1,intervals[2:(nintervals-1)],
                         intervals[nintervals]+.1))))
   }
   max.freq_max(sapply(freq,max))
   for (k in 1:nvectors)
   {
    line[[k]]_matrix(0,(nintervals-1),2)
    line[[k]][,1]_intervals[1:(nintervals-1)]
    line[[k]][,2]_intervals[2:nintervals]
    qlines[[k]]_c(line[[k]][1,1]-(line[[k]][1,2]-line[[k]][1,1])/2,
             line[[k]][(nintervals-1),2]+(line[[k]][1,2]-line[[k]][1,1])/2)
    midpoints[[k]]_line[[k]][,1]+(line[[k]][,2]-line[[k]][,1])/2
    upperedge[[k]]_low[k]+(freq[[k]]/max.freq)*boxwidth[k]
    loweredge[[k]]_rep(low[k],2)
    xpoints[[k]]_c(line[[k]][1,1]-(line[[k]][1,2]-line[[k]][1,1])/2,
                   midpoints[[k]],
                   line[[k]][(nintervals-1),2]+(line[[k]][1,2]-line[[k]][1,1])/2)
    ypoints[[k]]_c(low[k],upperedge[[k]],low[k])
   }
  }
  else
  {
   for (k in 1:nvectors)
   {
    if (all(is.na(breaks)))
    intervals_seq(min(x[[k]]),max(x[[k]]),
                       by=(max(x[[k]])-min(x[[k]]))/nclass)
    else intervals_breaks
    nintervals_length(intervals)
    freq[[k]]_as.vector(table(cut(x[[k]],
                breaks=c(intervals[1]-.1,intervals[2:(nintervals-1)],
                         intervals[nintervals]+.1))))
    line[[k]]_matrix(0,(nintervals-1),2)
    line[[k]][,1]_intervals[1:(nintervals-1)]
    line[[k]][,2]_intervals[2:nintervals]
    qlines[[k]]_c(line[[k]][1,1]-(line[[k]][1,2]-line[[k]][1,1])/2,
             line[[k]][(nintervals-1),2]+(line[[k]][1,2]-line[[k]][1,1])/2)
    midpoints[[k]]_line[[k]][,1]+(line[[k]][,2]-line[[k]][,1])/2
    upperedge[[k]]_low[k]+(freq[[k]]/max(freq[[k]]))*boxwidth[k]
    loweredge[[k]]_rep(low[k],2)
    xpoints[[k]]_c(line[[k]][1,1]-(line[[k]][1,2]-line[[k]][1,1])/2,
                   midpoints[[k]],
                   line[[k]][(nintervals-1),2]+(line[[k]][1,2]-line[[k]][1,1])/2)
    ypoints[[k]]_c(low[k],upperedge[[k]],low[k])
   }
  }
#-------------------------Draw frame and polygon----------------------------
  if (all(is.na(xlim))==T)
   fake.x_c(min(sapply(qlines,min)),max(sapply(qlines,max)))
  else fake.x_xlim
  plot(fake.x,c(0,1),type="n",axes=F,ylab="",xlab=xlab)
  if (axes==T)
  {
   axis(1,xaxt=xaxt)
   box()
  }
  for (k in 1:nvectors)
  {
   polygon(xpoints[[k]],ypoints[[k]],density=box.density,angle=box.angle,
           col=box.col,lwd=line.width,border=F)
   lines(xpoints[[k]],ypoints[[k]],lwd=line.width,col=line.col)
   lines(c(xpoints[[k]][1],xpoints[[k]][(nintervals-1)+2]),loweredge[[k]],
         lwd=line.width,col=line.col)
   if (box.label.on==T)
    text(midpoints[[k]],upperedge[[k]]+box.label.above,
         freqtext[[k]],cex=box.label.cex)
  }
}
#---------------------------------------------------------------------------
check.input.blip_function(x,nvectors,box.p,box.p.bar,graph.type,graph.width,
point.pattern,line.p,line.sd,line.se)
{
 nboxs_length(box.p)
#------------------ Check input error and put out error message------------
 if (length(line.sd)>1)
  stop("Only one value is needed for line.sd option")
 if (length(line.se)>1)
  stop("Only one value is needed for line.se option")
 if (graph.type != "b" & graph.type != "c" )
   stop("Incorrect value for graph.type")
 if (graph.width !="f" & graph.width !="v")
   stop("Incorrect value for graph.width")
 if (point.pattern !="on-line" & point.pattern !="stacking" &
      point.pattern !="evenly-spaced" & point.pattern !="jittered" &
      point.pattern !="max-range" & point.pattern !="vertical-bar")
  stop("Incorrect value for point.pattern")
#--------------------------------------------- Check p option;
 if ((any(is.na(box.p))==T) & nboxs>1)
  stop("Incorrect value in box.p option")
 if (all(!is.na(box.p))==T)                         # p!=NA
 {
  if (all(box.p[box.p!=0]==sort(box.p[box.p!=0]))==F)
   stop("Percentiles are not in an increasing order")
  if (nboxs<2)
   stop("Minimum two percentiles are needed for each box")
  if (any(duplicated(box.p[box.p!=0]))==T)
   stop("Duplicated percentiles are not allowed")
  if (any(box.p>1) | any(box.p<0))
   stop("box.p has to be between 0 and 1")
  if (box.p[nboxs]==0)
   stop("The last element of box.p cannot be 0")
  if (all(box.p==0))
   stop("All elements of box.p being 0 is not allowed")
  true.p_box.p
  if (box.p[1]==0 & box.p[2]==0) true.p[1]_true.p[2]_NA
  if (box.p[(nboxs-1)]==0 & box.p[nboxs]==1 & nboxs>2)
   true.p[(nboxs-1)]_true.p[nboxs]_NA
  true.p_true.p[!is.na(true.p)]
  real.breakpt_(2:length(true.p))[true.p[2:length(true.p)]==0]
  nreal.breakpt_length(real.breakpt)
  box.bounds_matrix(0,(nreal.breakpt+1),2)
  box.bounds[1,1]_true.p[1]
  if (nreal.breakpt==0)
  {
   box.bounds[1,2]_true.p[length(true.p)]
   if (box.bounds[1,1]==box.bounds[1,2])
    stop("At least two values for each box")
  }
  else
  {
   box.bounds[1,2]_true.p[real.breakpt[1]-1]
   if (box.bounds[1,1]==box.bounds[1,2])
    stop("At least two values for each box")
   box.bounds[nreal.breakpt+1,1]_true.p[real.breakpt[nreal.breakpt]+1]
   box.bounds[nreal.breakpt+1,2]_true.p[length(true.p)]
   if (box.bounds[nreal.breakpt+1,1]==box.bounds[nreal.breakpt+1,2])
    stop("At least two values for each box")
   if (nreal.breakpt>1)
   {
    for (i in 1:(nreal.breakpt-1))
    {
     box.bounds[i+1,1]_true.p[real.breakpt[i]+1]
     box.bounds[i+1,2]_true.p[real.breakpt[i+1]-1]
     if (box.bounds[i+1,1]==box.bounds[i+1,2])
      stop("Minimun two values needed for each box")
    }
   }
  }
 }
#---------------------------------------------Check box.p.bar option;
 if (any(is.numeric(box.p.bar)))
  stop("No numeric value is allowed for box.p.bar option")

 if (any(!is.na(box.p)))
 {
  if (length(box.p.bar) != nboxs)
   stop("Incorrect value for box.p.bar option")
 }
#---------------------------------------------Check line.p option;
 if (any(!is.na(line.p)))
 {
  nlines_length(line.p)
  if (all(line.p[line.p!=0]==sort(line.p[line.p!=0]))==F)
   stop("Percentiles are not in an increasing order")
  if (any(duplicated(line.p[line.p!=0]))==T)
   stop("Duplicated percentiles are not allowed")
  line.breakpt_(2:nlines)[line.p[2:nlines]==0]
  nline.breakpt_length(line.breakpt)
  line.bounds_matrix(0,nline.breakpt+1,2)
  line.bounds[1,1]_line.p[1]
  if (nline.breakpt==0)
  {
   line.bounds[1,2]_line.p[nlines]
   if (line.bounds[1,1]==line.bounds[1,2])
    stop("At least two values for each line segment")
  }
  else
  {
   line.bounds[1,2]_line.p[line.breakpt[1]-1]
   if (line.bounds[1,1]==line.bounds[1,2])
    stop("At least two values for each line segment")
   line.bounds[nline.breakpt+1,1]_line.p[line.breakpt[nline.breakpt]+1]
   line.bounds[nline.breakpt+1,2]_line.p[nlines]
   if (line.bounds[nline.breakpt+1,1]==line.bounds[nline.breakpt+1,2])
    stop("At least two values for each line segment")
   if (nline.breakpt>1)
   {
    for (j in 1:(nline.breakpt-1))
    {
     line.bounds[j+1,1]_line.p[line.breakpt[j]+1]
     line.bounds[j+1,2]_line.p[line.breakpt[j+1]-1]
     if (line.bounds[j+1,1]==line.bounds[j+1,2])
      stop("At least two values for each line segment")
    }
   }
  }
  if (all(!is.na(box.p))==T)
  {
   for (i in 1:dim(line.bounds)[1])
   {
    for (j in 1:dim(box.bounds)[1])
    {
     if (line.bounds[i,1]>box.bounds[j,1] & line.bounds[i,2]<box.bounds[j,2])
      stop("Incorrect values for line.p option")
     if (line.bounds[i,1]<box.bounds[j,1] & line.bounds[i,2]>box.bounds[j,2])
      stop("Incorrect values for line.p option")
     if (line.bounds[i,1]>box.bounds[j,1] & line.bounds[i,1]<box.bounds[j,2])
      stop("Incorrect values for line.p option")
     if (line.bounds[i,2]>box.bounds[j,1] & line.bounds[i,2]< box.bounds[j,2])
      stop("Incorrect values for line.p option")
    }
   }
  }
 }
#----------------------------Exclude points inside box and line------------
 xpoints_x
 if (any(!is.na(box.p)) | any(!is.na(line.p)))
 {
  for (k in 1:nvectors)
  {
   if (any(!is.na(line.p)))
   {
    for (i in 1:(nline.breakpt+1))
    {
     xpoints[[k]][((xpoints[[k]]>=quantile(x[[k]],probs=line.bounds[i,1])) &
                  (xpoints[[k]]<=quantile(x[[k]],probs=line.bounds[i,2])))]_NA
    }
    xpoints[[k]]_xpoints[[k]][!is.na(xpoints[[k]])]
   }
   if (any(!is.na(box.p)))
   {
    for (i in 1:(nreal.breakpt+1))
    {
     xpoints[[k]][((xpoints[[k]]>=quantile(x[[k]],probs=box.bounds[i,1])) &
                  (xpoints[[k]]<=quantile(x[[k]],probs=box.bounds[i,2])))]_NA
    }
    xpoints[[k]]_xpoints[[k]][!is.na(xpoints[[k]])]
   }
  }
 }
 if (all(is.na(box.p))) box.bounds_NA
  return(xpoints=xpoints,box.bounds=box.bounds)
}
#---------------------------------------------------------------------------
frame.x_function(x,nvectors,n,mean.pch,mean.cex,mean.col,line.sd,line.se,xlim,
                 graph.type,boxcenter,low,xpoints,line.width,
                 line.col,axes,xaxt,xlab)
{
#-----------------------------------------------------------------------------
#            Calculate sd, se and mean
#----------------------------------------------------------------------------
   if (!is.na(line.se))
   {
    mean.x_vector()
    se.x_vector()
    qline.se_list()
    for (k in 1:nvectors)
    {
     n[k]_length(x[[k]])
     mean.x[k]_mean(x[[k]])
     se.x[k]_sqrt(sum((x[[k]]-mean.x[k])^2)/(n[k]*(n[k]-1)))
     qline.se[[k]]_c(mean.x[k]-line.se*se.x[k],mean.x[k]+line.se*se.x[k])
     xpoints[[k]][(xpoints[[k]]>=qline.se[[k]][1]
                 & xpoints[[k]]<=qline.se[[k]][2])]_NA
     xpoints[[k]]_xpoints[[k]][!is.na(xpoints[[k]])]

    }
   }
   if (!is.na(line.sd))
   {
    mean.x_vector()
    sd.x_vector()
    qline.sd_list()
    for (k in 1:nvectors)
    {
     n[k]_length(x[[k]])
     mean.x[k]_mean(x[[k]])
     sd.x[k]_sqrt(var(x[[k]]))
     qline.sd[[k]]_c(mean.x[k]-line.sd*sd.x[k],mean.x[k]+line.sd*sd.x[k])
     xpoints[[k]][(xpoints[[k]]>=qline.sd[[k]][1]
                & xpoints[[k]]<=qline.sd[[k]][2])]_NA
     xpoints[[k]]_xpoints[[k]][!is.na(xpoints[[k]])]
    }
   }
   if (!is.na(mean.pch))
   {
    mean.x_vector()
    for (k in 1:nvectors)
    {
     mean.x[k]_mean(x[[k]])
    }
   }
#-----------------------------------------------------------------------------
#                  Build up basic frame
#---------------------------------------------------------------------------
  if (!is.na(line.sd) & !is.na(line.se))
  {
   if (all(is.na(xlim))==T)
   {
    fake1.x_c(min(sapply(qline.sd,min)),max(sapply(qline.sd,max)))
    fake2.x_c(min(sapply(qline.se,min)),max(sapply(qline.se,max)))
    x.min_min(sapply(x,min))
    x.max_max(sapply(x,max))
    if (x.min>=0) fake.x.min_.95*x.min
    else fake.x.min_1.05*x.min
    if (x.max>=0) fake.x.max_1.05*x.max
    else fake.x.max_.95*x.max
    fake3.x_c(fake.x.min,fake.x.max)
    fake.x_c(min(fake1.x[1],fake2.x[1],fake3.x[1]),
             max(fake1.x[2],fake2.x[2],fake3.x[2]))
   }
   else fake.x_xlim
  }
  else
  {
   if (!is.na(line.se))
   {
    if (all(is.na(xlim))==T)
    {
     fake1.x_c(min(sapply(qline.se,min)),max(sapply(qline.se,max)))
     x.min_min(sapply(x,min))
     x.max_max(sapply(x,max))
     if (x.min>=0) fake.x.min_.95*x.min
     else fake.x.min_1.05*x.min
     if (x.max>=0) fake.x.max_1.05*x.max
     else fake.x.max_.95*x.max
     fake2.x_c(fake.x.min,fake.x.max)
     fake.x_c(min(fake1.x[1],fake2.x[1]),max(fake1.x[2],fake2.x[2]))
    }
    else fake.x_xlim
   }
   else if (!is.na(line.sd))
   {
    if (all(is.na(xlim))==T)
    {
     fake1.x_c(min(sapply(qline.sd,min)),max(sapply(qline.sd,max)))
     x.min_min(sapply(x,min))
     x.max_max(sapply(x,max))
     if (x.min>=0) fake.x.min_.95*x.min
     else fake.x.min_1.05*x.min
     if (x.max>=0) fake.x.max_1.05*x.max
     else fake.x.max_.95*x.max
     fake2.x_c(fake.x.min,fake.x.max)
     fake.x_c(min(fake1.x[1],fake2.x[1]),max(fake1.x[2],fake2.x[2]))
    }
    else fake.x_xlim
   }
   else
   {
    if (all(is.na(xlim))==T)
    {
     x.min_min(sapply(x,min))
     x.max_max(sapply(x,max))
     if (x.min>=0) fake.x.min_.95*x.min
     else fake.x.min_1.05*x.min
     if (x.max>=0) fake.x.max_1.05*x.max
     else fake.x.max_.95*x.max
     fake.x_c(fake.x.min,fake.x.max)
    }
    else fake.x_xlim
   }
  }
  plot(fake.x,c(0,1),type="n",axes=F,ylab="",xlab=xlab)
  if (axes==T)
  {
   axis(1,xaxt=xaxt)
   box()
  }
#----------------------------Draw sd,se,mean-------------------------------
  if (!is.na(line.sd))
  {
   for (k in 1:nvectors)
   {
    if (graph.type=="c") lines(qline.sd[[k]],rep(boxcenter[k],2),
                               lwd=line.width,col=line.col)
    else lines(qline.sd[[k]],rep(low[k],2),lwd=line.width,col=line.col)
   }
  }
  if (!is.na(line.se))
  {
   for (k in 1:nvectors)
   {
    if (graph.type=="c") lines(qline.se[[k]],rep(boxcenter[k],2),
                               lwd=line.width,col=line.col)
    else lines(qline.se[[k]],rep(low[k],2),lwd=line.width,col=line.col)
   }
  }
  if (!is.na(mean.pch))
  {
   for (k in 1:nvectors)
   {
    if (graph.type=="c") points(mean.x[k],boxcenter[k],pch=mean.pch,
                                col=mean.col,cex=mean.cex)
    else points(mean.x[k],low[k],pch=mean.pch,col=mean.col,cex=mean.cex)
   }
  }
 return(xpoints=xpoints)
}
#----------------------------------------------------------------------------
draw.box_function(x,box.p,box.p.bar,box.bounds,nvectors,n,
                  graph.width,graph.type,graph.uniform,
                  point.pattern,nclass,boxcenter,boxwidth,
                  low,up,box.label.on,box.label.cex,box.label.above,
                  box.density,box.angle,box.col,line.width,line.col)
{
#-------------------------------------------------------------------------
# Determine the range of percentile at each point for drawing of box plot
#-------------------------------------------------------------------------
 nboxs_length(box.p)
 xunique_list()
 range.p.xunique_list()
 p.x_list()
 p.max_vector()
 delta_vector()
 if (graph.uniform==T)
 {
  x.max_max(sapply(x,max))
  x.min_min(sapply(x,min))
 }
 if (graph.width=="v")
 {
  for (k in 1:nvectors)
  {
   if (graph.uniform==T) delta[k]_(x.max-x.min)/(2*nclass)
   else  delta[k]_(max(x[[k]])-min(x[[k]]))/(2*nclass)
   xunique[[k]]_unique(x[[k]])
   xunique.l_xunique[[k]]-delta[k]
   xunique.u_xunique[[k]]+delta[k]
   ranks.x_sort.list(x[[k]])
   for (i in unique(x[[k]][duplicated(x[[k]])]))
   {
    ranks.x[x[[k]]==i]_max(ranks.x[x[[k]]==i])
   }
   p.x[[k]]_ranks.x/n[k]
   p.xunique.l_approx(x[[k]],p.x[[k]],xout=xunique.l)$y
   p.xunique.l[is.na(p.xunique.l)]_0
   p.xunique.u_approx(x[[k]],p.x[[k]],xout=xunique.u)$y
   p.xunique.u[is.na(p.xunique.u)]_1
   range.p.xunique[[k]]_p.xunique.u-p.xunique.l
   p.max[k]_max(range.p.xunique[[k]])
  }
  if (graph.uniform==T) p.max_rep(max(p.max),length(p.max))
 }
#----------------------------------------------------------------------------
#               Determine the height of the box
#----------------------------------------------------------------------------
 if (all(!is.na(box.p)==T))
 {
  for (k in 1:nvectors)
  {
   q_quantile(x[[k]],probs=box.p)
   nq_length(q)
   if (graph.width=="f")                        # Fixed width box
   {
    loweredge_rep(low[k],nq)
    upperedge_rep(up[k],nq)
   }
   else                                          # Variable width box
   {
    q.l_q-delta[k]
    q.u_q+delta[k]
    p.q.l_approx(x[[k]],p.x[[k]],xout=q.l)$y
    p.q.l[is.na(p.q.l)]_0
    p.q.u_approx(x[[k]],p.x[[k]],xout=q.u)$y
    p.q.u[is.na(p.q.u)]_1
    range.q_p.q.u-p.q.l
    range.q[range.q>p.max[k]]_p.max[k]
    if (graph.type=="c")                        # centered box
    {
     if (graph.uniform==T)
     {
      loweredge_boxcenter[k]-((range.q/p.max[k])*(n[k]/max(n))*boxwidth[k])/2
      upperedge_boxcenter[k]+((range.q/p.max[k])*(n[k]/max(n))*boxwidth[k])/2
     }
     else
     {
      loweredge_boxcenter[k]-((range.q/p.max[k])*boxwidth[k])/2
      upperedge_boxcenter[k]+((range.q/p.max[k])*boxwidth[k])/2
     }
    }
    else                                        # based box
    {
     loweredge_rep(low[k],nq)
     if (graph.uniform==T)
      upperedge_low[k]+(range.q/p.max[k])*(n[k]/max(n))*boxwidth[k]
     else
      upperedge_low[k]+(range.q/p.max[k])*boxwidth[k]
    }
   }
#------------------------------------------------------------------------------
#                  Draw the box
#-----------------------------------------------------------------------------
  for (i in 1:(dim(box.bounds)[1]))
  {
     q.box.bounds_quantile(x[[k]],probs=box.bounds[i,])
     q.box.bounds_q[q>=q.box.bounds[1] & q<=q.box.bounds[2]]
     hi_upperedge[box.p>=box.bounds[i,1] & box.p<=box.bounds[i,2]]
     lo_loweredge[box.p>=box.bounds[i,1] & box.p<=box.bounds[i,2]]
     polygon(c(q.box.bounds,rev(q.box.bounds)),
              c(hi,rev(lo)),border=F,
              density=box.density,angle=box.angle,col=box.col,lwd=line.width)
  }
  q[(2:nboxs)[box.p[2:nboxs]==0]]_NA
  qnames_names(q)
  if (all(is.na(box.p.bar))==T)
   segments(q,loweredge,q,upperedge,lwd=line.width,col=line.col)
  else if (any(box.p.bar)==T)
  {
   segments(q[box.p.bar==T],loweredge[box.p.bar==T],
            q[box.p.bar==T],upperedge[box.p.bar==T],lwd=line.width,
            col=line.col)
  }
  lines(q,loweredge,lwd=line.width,col=line.col)
  lines(q,upperedge,lwd=line.width,col=line.col)
  if (box.label.on==T)
   text(q[box.p.bar==T],upperedge[box.p.bar==T]+box.label.above,
        qnames[box.p.bar==T],cex=box.label.cex)
 }
}
return(xunique=xunique,range.p.xunique=range.p.xunique,p.max=p.max)
}
#-------------------------------------------------------------------------
draw.line_function(x,nvectors,line.p,graph.type,boxcenter,low,line.width,
                   line.col)
{
#----------------------------------------------------------------------------
#               Draw lines in line.p option
#----------------------------------------------------------------------------
 if (any(!is.na(line.p)))
 {
  qlines_list()
  nlines_length(line.p)
  for (k in 1:nvectors)
  {
   qlines[[k]]_quantile(x[[k]],probs=line.p)
   nqlines_length(qlines[[k]])
   qlines[[k]][(2:nlines)[line.p[2:nlines]==0]]_NA
   if (graph.type=="c") lines(qlines[[k]],rep(boxcenter[k],nqlines),
                              lwd=line.width,col=line.col)
   else  lines(qlines[[k]],rep(low[k],nqlines),lwd=line.width,col=line.col)
  }
 }
}
#----------------------------------------------------------------------------
draw.point_function(xpoints,nvectors,n,graph.width,graph.type,graph.uniform,
                    point.pattern,ans2,boxwidth,boxcenter,low,up,
                    seed,line.width,line.col,point.type,point.cex,
                    point.pch,point.col)
{
#--------------------------------------------------------------------------
#    Determine the value of y-axis of points that are to be drawn
#--------------------------------------------------------------------------
 if(any(!is.na(xpoints)))
 {
  if (point.pattern=="stacking")
  {
   tie.matrix_list()
   max.tie_vector()
   for (k in 1:nvectors)
   {
    xpoints.unique_unique(xpoints[[k]])
    tie.matrix[[k]]_matrix(NA,length(xpoints.unique),2)
    for (i in (1:length(xpoints.unique)))
    {
     tie.matrix[[k]][i,]_
     c(xpoints.unique[i],sum(xpoints[[k]]==xpoints.unique[i]))
    }
    max.tie[k]_max(tie.matrix[[k]][,2])
   }
   if (graph.uniform==T) increment_boxwidth/max(max.tie)
   else increment_boxwidth/max.tie
  }
  range.p.xunique_ans2$range.p.xunique
  p.max_ans2$p.max
  xunique_ans2$xunique
  for (k in 1:nvectors)
  {
   nxpoints_length(xpoints[[k]])
   if (point.pattern=="on-line")
   {
    if (graph.type=="b") ypoints_rep(low[k],nxpoints)
    else ypoints_rep(boxcenter[k],nxpoints)
   }
   else if (point.pattern=="jittered")
   {
    set.seed(seed)
    if (graph.width=="f") ypoints_runif(nxpoints,low[k],up[k])
    else
    {
     ypoints_vector()
     height.xpoints_vector()
     for (i in 1:nxpoints)
     {
      if (graph.uniform==T)
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==xpoints[[k]][i]]/p.max[k])*
       (n[k]/max(n))*boxwidth[k]
      else
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==xpoints[[k]][i]]/p.max[k])*
        boxwidth[k]
      if (graph.type=="b")
       ypoints[i]_runif(1,low[k],(height.xpoints[i]+low[k]))
      else ypoints[i]_runif(1,(boxcenter[k]-height.xpoints[i]/2),
                              (boxcenter[k]+height.xpoints[i]/2))
     }
    }
   }
   else if (point.pattern=="max-range")
   {
    ypoints_vector()
    height.xpoints_vector()
    for (i in 1:nxpoints)
    {
     if (graph.width=="f") ypoints[i]_up[k]
     else
     {
      if (graph.uniform==T)
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==xpoints[[k]][i]]/p.max[k])*
        (n[k]/max(n))*boxwidth[k]
      else
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==xpoints[[k]][i]]/p.max[k])*
        boxwidth[k]
      if (graph.type=="b")
       ypoints[i]_low[k]+height.xpoints[i]
      else ypoints[i]_boxcenter[k]+height.xpoints[i]/2
     }
    }
   }
   else if (point.pattern=="vertical-bar")
   {
    ypoints.up_vector()
    ypoints.low_vector()
    height.xpoints_vector()
    for (i in 1:nxpoints)
    {
     if (graph.width=="f")
     {
      ypoints.up[i]_up[k]
      ypoints.low[i]_low[k]
     }
     else
     {
      if (graph.uniform==T)
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==xpoints[[k]][i]]/p.max[k])*
        (n[k]/max(n))*boxwidth[k]
      else
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==xpoints[[k]][i]]/p.max[k])*
        boxwidth[k]
      if (graph.type=="b")
      {
       ypoints.up[i]_low[k]+height.xpoints[i]
       ypoints.low[i]_low[k]
      }
      else
      {
       ypoints.up[i]_boxcenter[k]+height.xpoints[i]/2
       ypoints.low[i]_boxcenter[k]-height.xpoints[i]/2
      }
     }
    }
   }
   else if (point.pattern=="evenly-spaced")
   {
    ypoints_vector()
    xpoints.unique_unique(xpoints[[k]])
    tie.matrix_matrix(NA,length(xpoints.unique),2)
    if (graph.width=="v") height.xpoints_vector()
    j_1
    for (i in (1:length(xpoints.unique)))
    {
     tie.matrix[i,]_c(xpoints.unique[i],sum(xpoints[[k]]==xpoints.unique[i]))
     if (graph.width=="v")
     {
      if (graph.uniform==T)
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==tie.matrix[i,1]]/p.max[k])*
       (n[k]/max(n))*boxwidth[k]
      else
       height.xpoints[i]_
       (range.p.xunique[[k]][xunique[[k]]==tie.matrix[i,1]]/p.max[k])*
        boxwidth[k]
     }
     if (tie.matrix[i,2]==1)
     {
      if (graph.width=="f")
      {
       if (graph.type=="c") ypoints[j]_boxcenter[k]
       else ypoints[j]_low[k]
      }
      else
      {
      if (graph.type=="c") ypoints[j]_boxcenter[k]
      else ypoints[j]_low[k]+height.xpoints[i]/2
      }
      j_j+1
     }
     else if (tie.matrix[i,2]>1)
     {
      for (l in (0:(tie.matrix[i,2]-1)))
      {
       if (graph.width=="f") ypoints[j]_low[k]+(boxwidth[k]/
                                        (tie.matrix[i,2]-1))*l
       else
       {
        if (graph.type=="c")
         ypoints[j]_(boxcenter[k]-height.xpoints[i]/2+(height.xpoints[i]/
                     (tie.matrix[i,2]-1))*l)
        else ypoints[j]_low[k]+(height.xpoints[i]/(tie.matrix[i,2]-1))*l
       }
       j_j+1
      }
     }
    }
   }
   else                                         # stacking
   {
    ypoints_vector()
    xpoints.unique_unique(xpoints[[k]])
    j_1
    for (i in (1:length(xpoints.unique)))
    {
     if (tie.matrix[[k]][i,2]==1)
     {
      if (graph.type=="b") ypoints[j]_low[k]
      else ypoints[j]_boxcenter[k]
      j_j+1
     }
     if (tie.matrix[[k]][i,2]>1)
     {
      for (l in (1:tie.matrix[[k]][i,2]))
      {
       if (graph.type=="b") ypoints[j]_low[k]+increment[k]*(l-1)
       else ypoints[j]_
       boxcenter[k]-increment[k]*(tie.matrix[[k]][i,2]-1)/2+increment[k]*(l-1)
       j_j+1
      }
     }
    }
  }
#----------------------------------------------------------------------------
#                  Draw the points
#----------------------------------------------------------------------------
  if (point.pattern!="vertical-bar")
   points(xpoints[[k]],ypoints,pch=point.pch,cex=point.cex,type=point.type,
          col=point.col)
  else
    segments(xpoints[[k]],ypoints.up,xpoints[[k]],ypoints.low,lwd=line.width,
             col=line.col)
 }#end of k loop
}
}