##############################################################################
##############################################################################
###   my_ruv_starter_analysis
##############################################################################
##############################################################################

my_ruv_starter_analysis <-
function(Y, X, ctl, Z=1, v=NULL, pctl=NULL, genecoloring=NULL, samplecoloring=NULL, genetexts=NULL, sampletexts=NULL, genesymbols=NULL, samplesymbols=NULL, geneinfo=NULL, rankbybeta=FALSE, topN=40, topcount_thresholds=c(20,40,60,80,100), rankset=NULL, kset=c(1,2,3,5,7,10,15,20,30,50,75,100,200,500,1000,2000,5000,10000), factorset=1:5, bin=10, do_general=TRUE, do_unadjusted=TRUE, do_ruv2=TRUE, do_ruv4=TRUE, do_ruvinv=TRUE, do_ruvrinv=TRUE, do_pptable=TRUE, outdir="html", initialize_collapsed=FALSE, webtitle="RUV Starter Analysis", inputcheck=TRUE, verbose=FALSE)
{
  # Sanity check on input variables
  if (verbose) cat("Input Check\n")
  if (inputcheck) inputcheck1(Y,X,Z,ctl)
  # Adjust for v, if given
  if (verbose) cat("adjust_v\n")
  Y = adjust_v(Y,v,ctl)
  ## Here we set up some variables.  
  ## We define m, n, nc, Yc, Z, q, kset, fullctl, rankset, genecoloring, and W_ctl and W_fullctl.  
  ## W_ctl are the left singular values of Yc after projecting away X and Z.  
  ## W_fullctl are the left singular values of Y after projecting away X and Z.
  if (verbose) cat("Setup\n")
  m = nrow(Y)
  n = ncol(Y)
  nc = sum(ctl)
  Yc = Y[,ctl,drop=FALSE]
  if (ncol(X) > 1) 
  {
    cat("This routine is not meant for X with more than one column.\n")
    return()
  }
  if (is.null(Z))
  {
    q = 0
  }
  else if (length(Z) == 1) {if (Z==1)  
  {
    Z = matrix(1,m,1)
    q = 1
  }}
  else 
  {
    q = ncol(Z)
  }
  kset = kset[kset < m-q-1] # ignore values of K that are out of range
  factorset = factorset[factorset < m-q-1] # ignore factors that are out of range
  fullctl = rep(T,n)
  if (is.null(rankset)) rankset = rep(TRUE, n)
  topcount_thresholds = topcount_thresholds[topcount_thresholds < sum(rankset)] # ignore values that are out of range
  if (is.null(genecoloring))
  {
    genecoloring = rep(NA,n)
    genecoloring[ctl] = "green4"
    if (!is.null(pctl)) genecoloring[pctl] = "purple"
  }
  unadj = RUV2(Y, X, ctl, 0, Z=Z, inputcheck=FALSE)
  W_ctl = unadj$fullW 
  W_fullctl = RUV2(Y, X, fullctl, 0, Z=Z, inputcheck=FALSE)$fullW

  ## Here we read in some text files that are needed in the production of the web page (labels, etc.)
  if (verbose) cat("read misc files\n")
  filelist = scan(paste(system.file(package="ruv.extras"), "misc/filelist.txt", sep="/"), what="character", quiet=TRUE)
  combinedheaders = scan(paste(system.file(package="ruv.extras"), "misc/combinedheaders.txt", sep="/"), what="character", sep="\n", quiet=TRUE)
  filestokeep = c(1,2,33,3:32)
  filelist = filelist[filestokeep]
  combinedheaders = combinedheaders[filestokeep]

  ################################################################################################
  ################################################################################################
  
  # Now we begin the "real work"

  hstart(overwrite=TRUE, basedir=outdir, title=webtitle, bodyfiles=paste(system.file(package="ruv.htmllatex"), "misc/javascript.js", sep="/")) # open the web page for output

  ###########################################
  ## General Info
  ###########################################
  
  # This is the first section of the web page.  

  if (verbose) cat("General Info\n")
  if(do_general)
  {
    hsectionstart("generalsec", "General Info", initialize_collapsed=initialize_collapsed) #This defines the section header / section in the web page
    my_general_info(m, n, nc, Y, Yc, X, Z, ctl, W_fullctl, W_ctl, samplecoloring, genecoloring, sampletexts, genetexts, samplesymbols, genesymbols, initialize_collapsed=initialize_collapsed)
    hsectionstop()
  }

  ###########################################
  ## Unadjusted
  ###########################################

  if (verbose) cat("Unadjusted\n")
  if(do_unadjusted)
  {
    hsectionstart("unadjsec", "Unadjusted", initialize_collapsed=initialize_collapsed)
    fit = variance_adjust(unadj, bin=bin)
    my_individual_analysis(fit, "unadj", Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, do_rle=TRUE, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
    hsectionstop()
  }

  ###########################################
  ## RUV-2 Analyses
  ###########################################

  if (verbose) cat("RUV-2\n")
  if(do_ruv2)
  {
    ## Combined RUV-2 Analyses
    if (verbose) cat("Combined\n")
    hsectionstart("combruv2sec", "RUV-2 Combined Analyses", initialize_collapsed=initialize_collapsed)
    if (!is.null(pctl))  counttables = my_ruv2toprankcounttables(Y, X, ctl, pctl, Z=Z, bin=bin, thresholds=topcount_thresholds, rankset=rankset, inputcheck=FALSE)
    else counttables = NULL
    my_combined_analysis(counttables, filelist, combinedheaders, "ruv2", kset, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
    hsectionstop()

    ## Individual RUV-2 Analyses
    if (verbose) cat("Individual\n")
    hsectionstart("ruv2sec", "RUV-2 Individual Analyses", initialize_collapsed=initialize_collapsed)
    W = RUV2(Y,X,ctl,0,Z=Z,inputcheck=FALSE)$fullW
    for (k in kset)
    {
      if (verbose) cat(k, ", ", sep="")
      hcat("<br>\n")
      hsubsectionstart(paste("ruv2k", k, "sec", sep=""), paste("RUV-2 K =", k), initialize_collapsed=initialize_collapsed)
      hcat("<br>\n")
      fit = RUV2(Y,X,ctl,k,Z=Z,fullW=W,inputcheck=FALSE)
      fit = variance_adjust(fit, bin=bin)
      my_individual_analysis(fit, "ruv2", Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, k=k, do_rle=TRUE, do_projection=TRUE, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
      hsubsectionstop()
    }
    if (verbose) cat("\n")
    hsectionstop()
  }

  ###########################################
  ## RUV-4 Analyses
  ###########################################

  if (verbose) cat("RUV-4\n")
  if(do_ruv4)
  {
    ## Combined RUV-4 Analyses
    if (verbose) cat("Combined\n")
    hsectionstart("combruv4sec", "RUV-4 Combined Analyses", initialize_collapsed=initialize_collapsed)
    if (!is.null(pctl))  counttables = my_ruv4toprankcounttables(Y, X, ctl, pctl, Z=Z, bin=bin, thresholds=topcount_thresholds, rankset=rankset, inputcheck=FALSE)
    else counttables = NULL
    my_combined_analysis(counttables, filelist, combinedheaders, "ruv4", kset, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
    hsectionstop()
  
    ## Individual RUV-4 Analyses
    if (verbose) cat("Individual\n")
    hsectionstart("ruv4sec", "RUV-4 Individual Analyses", initialize_collapsed=initialize_collapsed)
    W0 = RUV4(Y,X,ctl,0,Z=Z,inputcheck=FALSE)$fullW0
    for (k in kset)
    {
      if (verbose) cat(k, ", ", sep="")
      hsubsectionstart(paste("ruv4k", k, "sec", sep=""), paste("RUV-4 K =", k), initialize_collapsed=initialize_collapsed)
      fit = RUV4(Y,X,ctl,k,Z=Z,fullW0=W0,inputcheck=FALSE)
      fit = variance_adjust(fit, bin=bin)
      my_individual_analysis(fit, "ruv4", Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, k=k, do_rle=TRUE, do_projection=TRUE, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
      hsubsectionstop()
    }
    if (verbose) cat("\n")
    hsectionstop()
  }

  ###########################################
  ## RUV-inv
  ###########################################

  if (verbose) cat("RUV-inv\n")
  if(do_ruvinv)
  {
    hsectionstart("ruvinvsec", "RUV-inv", initialize_collapsed=initialize_collapsed)
    fit = RUVinv(Y,X,ctl,Z=Z,inputcheck=FALSE)
    fit = variance_adjust(fit, bin=bin)
    my_individual_analysis(fit, "ruvinv", Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, do_projection=TRUE, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
    hsectionstop() 
  }

  ###########################################
  ## RUV-rinv
  ###########################################

  if (verbose) cat("RUV-rinv\n")
  if(do_ruvrinv)
  {
    hsectionstart("ruvrinvsec", "RUV-rinv", initialize_collapsed=initialize_collapsed)
    fit = RUVrinv(Y,X,ctl,Z=Z,inputcheck=FALSE)
    fit = variance_adjust(fit, bin=bin)
    my_individual_analysis(fit, "ruvrinv", Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, do_projection=TRUE, rankbybeta=rankbybeta, initialize_collapsed=initialize_collapsed)
    hsectionstop()
  }

  ###########################################
  ## Projection Plot Table
  ###########################################

  if (verbose) cat("Projection Plot Table\n")
  if(do_pptable)
  {
    hsectionstart("projplottablesec", "Projection Plot Table", initialize_collapsed=initialize_collapsed)
    fullfit = RUVinv(Y,X,ctl,Z=Z,inputcheck=FALSE)
    my_projplottable(fullfit, c(0,kset), factorset, genecoloring, genetexts, genesymbols, initialize_collapsed=initialize_collapsed)
    hsectionstop()
  }

  hstop() # Close the web page
}





######################################################################################
## General Info
######################################################################################

my_general_info = function(m, n, nc, Y, Yc, X, Z, ctl, W_fullctl, W_ctl, samplecoloring, genecoloring, sampletexts, genetexts, samplesymbols, genesymbols, initialize_collapsed=FALSE)
{
  hcat("<font size=5> m = ", m, ", n = ", n, ", n_c = ", nc, "</font><br><br>\n")

  exprs = filenames = matrix("",1,2)
  filenames[1,1] = "screeY.png"
  filenames[1,2] = "screeYc.png"
  exprs[1,1] = "my_ruv_scree(Y, Z=Z)"
  exprs[1,2] = "my_ruv_scree(Yc, Z=Z)"
  hdynamicplotmat("Scree Plots", "scree", exprs, filenames, headers=c("Y", "Y_c"), dir="general", initialize_collapsed=initialize_collapsed)    

  exprs = filenames = matrix("",1,2)
  filenames[1,1] = "rleY.png"
  filenames[1,2] = "rleYc.png"
  exprs[1,1] = "my_ruv_rle(Y, XZW=NULL, coloring=samplecoloring)"
  exprs[1,2] = "my_ruv_rle(Yc, XZW=NULL, coloring=samplecoloring)"
  hdynamicplotmat("RLE Plots", "rle", exprs, filenames, headers=c("Y", "Y_c"), dir="general", initialize_collapsed=initialize_collapsed)

  exprs = filenames = matrix("",1,1)
  filenames[1,1] = "cancor.png"
  exprs[1,1] = "my_ruv_cancorplot(X, W_fullctl, W2=W_ctl, col=c(1, genecoloring[ctl][1]))"
  hdynamicplotmat("Cancor Plots", "cancor", exprs, filenames, dir="general", initialize_collapsed=initialize_collapsed)

  exprs = filenames = matrix("",1,1)
  filenames[1,1] = "pcgridplot.png"
  exprs[1,1] = "my_ruv_pcgridplot(W_fullctl, W2=W_ctl, kset=c(1,2,3), coloring=samplecoloring, symbols=samplesymbols, texts=sampletexts)"
  hdynamicplotmat("PC Plots", "pcgrid", exprs, filenames, dir="general", height=1000, width=1000, initialize_collapsed=initialize_collapsed)

  exprs = filenames = matrix("",1,1)
  filenames[1,1] = "alphagridplot.png"
  exprs[1,1] = "my_ruv_pcgridplot(t(Y)%*%W_fullctl[,1:3], W2=t(Y)%*%W_ctl[,1:3], kset=c(1,2,3), coloring=genecoloring, texts=genetexts, symbols=genesymbols)"
  hdynamicplotmat("Alpha Plots", "alphagrid", exprs, filenames, dir="general", height=1000, width=1000, initialize_collapsed=initialize_collapsed)
}

######################################################################################
## Individual Analysis
######################################################################################

my_individual_analysis <-
function(fit, id, Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo=NULL, topN, rankset=rankset, k=NULL, do_rle=FALSE, do_pval=TRUE, do_ctlpval=TRUE, do_projection=FALSE, do_variance=TRUE, do_counttable=TRUE, do_ranktable=TRUE, rankbybeta=FALSE, initialize_collapsed=FALSE)
{
  if (is.null(k)) kstring = ""
  else kstring = paste(" (K = ", k, ")", sep="")
  if (is.null(k)) kstring2 = ""
  else kstring2 = paste("k", k, sep="")
  if (is.null(rankset)) rankset = rep(TRUE,ncol(Y))
  if (do_rle)
  {
      exprs = filenames = matrix("",1,2)
      filenames[1,1] = "rleY.png"
      filenames[1,2] = "rleYc.png"
      exprs[1,1] = "my_ruv_rle(Y, XZW=cbind(fit$X,fit$Z,fit$W), coloring=samplecoloring)"
      exprs[1,2] = "my_ruv_rle(Y[,fit$ctl,drop=FALSE], XZW=cbind(fit$X,fit$Z,fit$W), coloring=samplecoloring)"
      hdynamicplotmat(paste("RLE Plots", kstring, sep=""), paste(id, kstring2, "rle", sep=""), exprs, filenames, headers=c("Y", "Y_c"), dir=id, initialize_collapsed=initialize_collapsed)
  }
  if (do_pval)
  {
      exprs = filenames = matrix("",3,5)
      filenames[1,1] = "unadjhist.png"
      filenames[1,2] = "ebayeshist.png"
      filenames[1,3] = "rsvarhist.png"
      filenames[1,4] = "rsvarebayeshist.png"
      filenames[1,5] = "evarhist.png"
      filenames[2,1] = "unadjplot.png"
      filenames[2,2] = "ebayesplot.png"
      filenames[2,3] = "rsvarplot.png"
      filenames[2,4] = "rsvarebayesplot.png"
      filenames[2,5] = "evarplot.png"
      filenames[3,1] = "unadjlogplot.png"
      filenames[3,2] = "ebayeslogplot.png"
      filenames[3,3] = "rsvarlogplot.png"
      filenames[3,4] = "rsvarebayeslogplot.png"
      filenames[3,5] = "evarlogplot.png"
      exprs[1,1] = "my_ruv_pvalhist(fit$p)"
      exprs[1,2] = "my_ruv_pvalhist(fit$p.ebayes)"
      exprs[1,3] = "my_ruv_pvalhist(fit$p.rsvar)"
      exprs[1,4] = "my_ruv_pvalhist(fit$p.rsvar.ebayes)"
      exprs[1,5] = "my_ruv_pvalhist(fit$p.evar)"
      exprs[2,1] = "my_ruv_pvalplot(fit$p, coloring=genecoloring, symbols=genesymbols, texts=genetexts)"
      exprs[2,2] = "my_ruv_pvalplot(fit$p.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts)"
      exprs[2,3] = "my_ruv_pvalplot(fit$p.rsvar, coloring=genecoloring, symbols=genesymbols, texts=genetexts)"
      exprs[2,4] = "my_ruv_pvalplot(fit$p.rsvar.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts)"
      exprs[2,5] = "my_ruv_pvalplot(fit$p.evar, coloring=genecoloring, symbols=genesymbols, texts=genetexts)"
      exprs[3,1] = "my_ruv_pvalplot(fit$p, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)"
      exprs[3,2] = "my_ruv_pvalplot(fit$p.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)"
      exprs[3,3] = "my_ruv_pvalplot(fit$p.rsvar, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)"
      exprs[3,4] = "my_ruv_pvalplot(fit$p.rsvar.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)"
      exprs[3,5] = "my_ruv_pvalplot(fit$p.evar, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)"
      hdynamicplotmat(paste("P-values", kstring, sep=""), paste(id, kstring2, "pval", sep=""), exprs, filenames, headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar"), dir=id, initialize_collapsed=initialize_collapsed)
  }
  if (do_ctlpval)
  {
      exprs = filenames = matrix("",3,5)
      filenames[1,1] = "unadjhist.png"
      filenames[1,2] = "ebayeshist.png"
      filenames[1,3] = "rsvarhist.png"
      filenames[1,4] = "rsvarebayeshist.png"
      filenames[1,5] = "evarhist.png"
      filenames[2,1] = "unadjplot.png"
      filenames[2,2] = "ebayesplot.png"
      filenames[2,3] = "rsvarplot.png"
      filenames[2,4] = "rsvarebayesplot.png"
      filenames[2,5] = "evarplot.png"
      filenames[3,1] = "unadjlogplot.png"
      filenames[3,2] = "ebayeslogplot.png"
      filenames[3,3] = "rsvarlogplot.png"
      filenames[3,4] = "rsvarebayeslogplot.png"
      filenames[3,5] = "evarlogplot.png"
      exprs[1,1] = "my_ruv_pvalhist(fit$p[fit$ctl])"
      exprs[1,2] = "my_ruv_pvalhist(fit$p.ebayes[fit$ctl])"
      exprs[1,3] = "my_ruv_pvalhist(fit$p.rsvar[fit$ctl])"
      exprs[1,4] = "my_ruv_pvalhist(fit$p.rsvar.ebayes[fit$ctl])"
      exprs[1,5] = "my_ruv_pvalhist(fit$p.evar[fit$ctl])"
      exprs[2,1] = "my_ruv_pvalplot(fit$p[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])"
      exprs[2,2] = "my_ruv_pvalplot(fit$p.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])"    
      exprs[2,3] = "my_ruv_pvalplot(fit$p.rsvar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])"
      exprs[2,4] = "my_ruv_pvalplot(fit$p.rsvar.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])"
      exprs[2,5] = "my_ruv_pvalplot(fit$p.evar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])"
      exprs[3,1] = "my_ruv_pvalplot(fit$p[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)"
      exprs[3,2] = "my_ruv_pvalplot(fit$p.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)"
      exprs[3,3] = "my_ruv_pvalplot(fit$p.rsvar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)"
      exprs[3,4] = "my_ruv_pvalplot(fit$p.rsvar.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)"
      exprs[3,5] = "my_ruv_pvalplot(fit$p.evar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)"
      hdynamicplotmat(paste("P-values (Controls)", kstring, sep=""), paste(id, kstring2, "pvalctl", sep=""), exprs, filenames, headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar"), dir=id, initialize_collapsed=initialize_collapsed)
  }
  if (do_projection)
  {
      exprs = filenames = matrix("",1,6)
      filenames[1,1] = "projection.png"
      filenames[1,2] = "topprojection.png"
      filenames[1,3] = "topprojection.ebayes.png"
      filenames[1,4] = "topprojection.rsvar.png"
      filenames[1,5] = "topprojection.rsvar.ebayes.png"
      filenames[1,6] = "topprojection.evar.png"
      if (rankbybeta) filenames = cbind(filenames, "topprojection.beta.png")
      exprs[1,1] = "my_ruv_projectionplot(fit$byx, fit, fit$bwx, coloring=genecoloring, symbols=genesymbols, texts=genetexts)"
      exprs[1,2] = "my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p, rankset=rankset)"
      exprs[1,3] = "my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.ebayes, rankset=rankset)"
      exprs[1,4] = "my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.rsvar, rankset=rankset)"
      exprs[1,5] = "my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.rsvar.ebayes, rankset=rankset)"
      exprs[1,6] = "my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.evar, rankset=rankset)"
      if (rankbybeta) exprs = cbind(exprs, "my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, -abs(fit$betahat), rankset=rankset)")
      headers = c("Standard Coloring", "standard", "ebayes", "rsvar", "rsvar ebayes", "evar")
      if (rankbybeta) headers = c(headers, "betahat")
      hdynamicplotmat(paste("Projection Plots", kstring, sep=""), paste(id, kstring2, "projplots", sep=""), exprs, filenames, headers=headers, dir=id, initialize_collapsed=initialize_collapsed)
}
  if (do_variance)
  {
      exprs = filenames = matrix("",1,6)
      filenames[1,1] = "varianceroot.png"
      filenames[1,2] = "topvarianceroot.png"  
      filenames[1,3] = "topvariancerootebayes.png"
      filenames[1,4] = "topvariancerootrsvar.png"
      filenames[1,5] = "topvariancerootrsvarebayes.png"
      filenames[1,6] = "topvariancerootevar.png"
      if (rankbybeta) filenames = cbind(filenames, "topvariancerootbeta.png")
      exprs[1,1] = "my_ruv_varianceplot(fit, fit$varbetahat, coloring=genecoloring, symbols=genesymbols, texts=genetexts, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\")"
      exprs[1,2] = "my_ruv_topvarianceplot(fit, fit$p, fit$varbetahat, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\", rankset=rankset)"
      exprs[1,3] = "my_ruv_topvarianceplot(fit, fit$p.ebayes, fit$varbetahat.ebayes, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\", rankset=rankset)"
      exprs[1,4] = "my_ruv_topvarianceplot(fit, fit$p.rsvar, fit$varbetahat.rsvar, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\", rankset=rankset)"
      exprs[1,5] = "my_ruv_topvarianceplot(fit, fit$p.rsvar.ebayes, fit$varbetahat.rsvar.ebayes, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\", rankset=rankset)"
      exprs[1,6] = "my_ruv_topvarianceplot(fit, fit$p.evar, fit$varbetahat.evar, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\", rankset=rankset)"
      if (rankbybeta) exprs = cbind(exprs, "my_ruv_topvarianceplot(fit, -abs(fit$betahat), fit$varbetahat, trans=my_fourthroot, xlab=\"(sigma^2 hat)^(1/4)\", ylab=\"(betahat^2)^(1/4)\", rankset=rankset)")
      headers = c("Standard Coloring", "standard", "ebayes", "rsvar", "rsvar ebayes", "evar")
      if (rankbybeta) headers = c(headers, "betahat")
      hdynamicplotmat(paste("Variance Plots", kstring, sep=""), paste(id, kstring2, "varplots", sep=""), exprs, filenames, headers=headers, dir=id, initialize_collapsed=initialize_collapsed)
  }
  if (do_counttable)
  {
    if (!is.null(pctl))
    {
      rankcount = my_ruv_toprankcount(fit$p[1,rankset], pctl[rankset], topcount_thresholds)
      rankcount.ebayes = my_ruv_toprankcount(fit$p.ebayes[1,rankset], pctl[rankset], topcount_thresholds)
      rankcount.rsvar = my_ruv_toprankcount(fit$p.rsvar[1,rankset], pctl[rankset], topcount_thresholds)
      rankcount.rsvar.ebayes = my_ruv_toprankcount(fit$p.rsvar.ebayes[1,rankset], pctl[rankset], topcount_thresholds)
      rankcount.evar = my_ruv_toprankcount(fit$p.evar[1,rankset], pctl[rankset], topcount_thresholds)
      rankcount.beta = my_ruv_toprankcount(-abs(fit$betahat[1,rankset]), pctl[rankset], topcount_thresholds)
      if (rankbybeta) hdynamictablerownest("Number of Top-ranked Positive Controls", paste(id, kstring2, "counts", sep=""), list(rankcount, rankcount.ebayes, rankcount.rsvar, rankcount.rsvar.ebayes, rankcount.evar, rankcount.beta), headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar", "betahat"), google=FALSE, initialize_collapsed=initialize_collapsed)
      else hdynamictablerownest("Number of Top-ranked Positive Controls", paste(id, kstring2, "counts", sep=""), list(rankcount, rankcount.ebayes, rankcount.rsvar, rankcount.rsvar.ebayes, rankcount.evar), headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar"), google=FALSE, initialize_collapsed=initialize_collapsed)
    }
  }
  if (do_ranktable)
  {
      geneinfo = cbind(geneinfo, index = 1:length(fit$p))
      ranktable = my_ruv_topranktable(fit$p, beta=fit$betahat, info=geneinfo, N=topN, rankset=rankset)
      ranktable.ebayes = my_ruv_topranktable(fit$p.ebayes, beta=fit$betahat, info=geneinfo, N=topN, rankset=rankset)
      ranktable.rsvar = my_ruv_topranktable(fit$p.rsvar, beta=fit$betahat, info=geneinfo, N=topN, rankset=rankset)
      ranktable.rsvar.ebayes = my_ruv_topranktable(fit$p.rsvar.ebayes, beta=fit$betahat, info=geneinfo, N=topN, rankset=rankset)
      ranktable.evar = my_ruv_topranktable(fit$p.evar, beta=fit$betahat, info=geneinfo, N=topN, rankset=rankset)
      ranktable.beta = my_ruv_topranktable(fit$p, beta=fit$betahat, rankby=-abs(fit$betahat), info=geneinfo, N=topN, rankset=rankset)
      ranktable.ebayes.beta = my_ruv_topranktable(fit$p.ebayes, beta=fit$betahat, rankby=-abs(fit$betahat), info=geneinfo, N=topN, rankset=rankset)
      ranktable.rsvar.beta = my_ruv_topranktable(fit$p.rsvar, beta=fit$betahat, rankby=-abs(fit$betahat), info=geneinfo, N=topN, rankset=rankset)
      ranktable.rsvar.ebayes.beta = my_ruv_topranktable(fit$p.rsvar.ebayes, beta=fit$betahat, rankby=-abs(fit$betahat), info=geneinfo, N=topN, rankset=rankset)
      ranktable.evar.beta = my_ruv_topranktable(fit$p.evar, beta=fit$betahat, rankby=-abs(fit$betahat), info=geneinfo, N=topN, rankset=rankset)
      hdynamictablerownest(paste("Top-ranked Genes", kstring, sep=""), paste(id, kstring2, "ranks", sep=""), list(ranktable, ranktable.ebayes, ranktable.rsvar, ranktable.rsvar.ebayes, ranktable.evar), headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar"), google=TRUE, initialize_collapsed=initialize_collapsed)
      if (rankbybeta) hdynamictablerownest(paste("Top-ranked Genes, Ranked by betahat", kstring, sep=""), paste(id, kstring2, "betaranks", sep=""), list(ranktable.beta, ranktable.ebayes.beta, ranktable.rsvar.beta, ranktable.rsvar.ebayes.beta, ranktable.evar.beta), headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar"), google=TRUE, initialize_collapsed=initialize_collapsed)
  }
}

######################################################################################
## Combined Analysis
######################################################################################

my_combined_analysis <-
function(counttables, filelist, combinedheaders, id, kset, rankbybeta=FALSE, initialize_collapsed=FALSE)
{
    if (!is.null(counttables))
    {
      exprs = filenames = matrix("",1,5)
      filenames[1,1] = "countplot.png"
      filenames[1,2] = "countplot.ebayes.png"
      filenames[1,3] = "countplot.rsvar.png"
      filenames[1,4] = "countplot.rsvar.ebayes.png"
      filenames[1,5] = "countplot.evar.png"
      if (rankbybeta) filenames = cbind(filenames, "countplot.beta.png")
      exprs[1,1] = "my_ruv_toprankcountplot(counttables$countmat, counttables$kset, counttables$thresholds)"
      exprs[1,2] = "my_ruv_toprankcountplot(counttables$countmat.ebayes, counttables$kset, counttables$thresholds)"
      exprs[1,3] = "my_ruv_toprankcountplot(counttables$countmat.rsvar, counttables$kset, counttables$thresholds)"
      exprs[1,4] = "my_ruv_toprankcountplot(counttables$countmat.rsvar.ebayes, counttables$kset, counttables$thresholds)"
      exprs[1,5] = "my_ruv_toprankcountplot(counttables$countmat.evar, counttables$kset, counttables$thresholds)"
      if (rankbybeta) exprs = cbind(exprs, "my_ruv_toprankcountplot(counttables$countmat.beta, counttables$kset, counttables$thresholds)")
      headers=c("standard", "ebayes", "rsvar", "rsvar ebayes", "evar")
      if (rankbybeta) headers = c(headers, "betahat")
      hdynamicplotmat("Top-ranked Count Plots", paste(id, "krankcount", sep=""), exprs, filenames, headers=headers, dir=id, initialize_collapsed=initialize_collapsed)
    }

    filemat = matrix("", length(kset), length(filelist))
    for (i in 1:nrow(filemat)) for (j in 1:ncol(filemat))
    {
      filemat[i,j] = paste("<img src=\"", id, "/", id, "k", kset[i], filelist[j], "\" alt=\"image\" />", sep="")
    }
    filemat = cbind(kset,filemat)
    hdynamictablemat("Combined Diagnostic Plots", paste(id, "combined", sep=""), filemat, headers = c("K",combinedheaders), initialize_collapsed=initialize_collapsed)
}

######################################################################################
## Projection Plot Table
######################################################################################

my_projplottable <- function(fullfit, kset, factorset, genecoloring, genetexts, genesymbols, initialize_collapsed=initialize_collapsed)
{
  exprs = filenames = matrix("", length(kset), length(factorset))
  for (i in 1:length(kset)) for (j in 1:length(factorset))
  {
    filenames[i,j] = paste("factorpp_k", kset[i], "f", factorset[j], ".png", sep="")
    exprs[i,j] = paste("my_ruv_factor_projectionplot(fullfit, ", kset[i], ", ", factorset[j], ", coloring=genecoloring, symbols=genesymbols, texts=genetexts)", sep="")
  }
  hdynamicplotmat("Factor Projection Plots", "fpp", exprs, filenames, headers=c(factorset), dir="proj_plots", initialize_collapsed=initialize_collapsed)
}