##############################################################################
##############################################################################
###   my_ruv_starter_analysis (simplest)
##############################################################################
##############################################################################

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, 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)
{
  # Adjust for v, if given
  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.
  m = nrow(Y)	# number of samples
  n = ncol(Y)	# number of features (genes)
  nc = sum(ctl)	# number of control features
  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)) # if no genecoloring is specified, the default is to set negative controls green and positive controls purple
  {
    genecoloring = rep(NA,n)
    genecoloring[ctl] = "green4"
    if (!is.null(pctl)) genecoloring[pctl] = "purple"
  }
  unadj = RUV2(Y, X, ctl, 0, Z=Z) # setting k=0 is equivalent to an unadjusted analysis (no RUV)
  W_ctl = unadj$fullW # W_ctl is the left singular vectors of Y_c
  W_fullctl = RUV2(Y, X, fullctl, 0, Z=Z)$fullW # W_fullctl is the left singular vectors of Y

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

  ###########################################
  ## General Info
  ###########################################
  
  # this generates the plots found in the "General Info" section of the web page
  if (do_general)
  {
    my_general_info(m, n, nc, Y, Yc, X, Z, ctl, W_fullctl, W_ctl, samplecoloring, genecoloring, sampletexts, genetexts, samplesymbols, genesymbols) 
  }

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

  if(do_unadjusted)
  {
    # calculate ebayes, rescaled, and empirical variances, along with associated p-values
    fit = variance_adjust(unadj, bin=bin) 
    # this generates the plots found in the "Unadjusted" section of the web page.  RLE plots are not generated by default, since they are not always appropriate (e.g. RUV-inv).  Therefore we must explicitly set do_rle=TRUE.
    my_individual_analysis(fit, Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, do_rle=TRUE) 

  }

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

  if(do_ruv2)
  {
    ## Combined RUV-2 Analyses
    # if there are positive controls available, count how many of them are top-ranked for various values of K.  "topcount_thresholds" specifies which values of K to use.
    if (!is.null(pctl))  counttables = my_ruv2toprankcounttables(Y, X, ctl, pctl, Z=Z, bin=bin, thresholds=topcount_thresholds, rankset=rankset) 
    else counttables = NULL
    # This produces the plots.
    my_combined_analysis(counttables)

    ## Individual RUV-2 Analyses
    W = RUV2(Y,X,ctl,0,Z=Z)$fullW  # only calculate W once, to speed things up
    for (k in kset)
    {
      # run the RUV2 analysis
      fit = RUV2(Y,X,ctl,k,Z=Z,fullW=W) 
      # calculate ebayes, rescaled, and empirical variances, along with associated p-values
      fit = variance_adjust(fit, bin=bin) 
      # this generates the plots found in the "RUV-2" section of the web page.  Projection plots are not generated by default, since they are not always appropriate (e.g. unadjusted analysis).  Therefore we must explicitly set do_projection=TRUE.  As above, we also must set do_RLE=TRUE.
      my_individual_analysis(fit, Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, k=k, do_rle=TRUE, do_projection=TRUE)
    }
  }

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

  # Analogous to RUV-2 section.  See comments above.
  if(do_ruv4) 
  {
    if (!is.null(pctl))  counttables = my_ruv4toprankcounttables(Y, X, ctl, pctl, Z=Z, bin=bin, thresholds=topcount_thresholds, rankset=rankset)
    else counttables = NULL
    my_combined_analysis(counttables)
  
    W0 = RUV4(Y,X,ctl,0,Z=Z)$fullW0 
    for (k in kset)
    {
      fit = RUV4(Y,X,ctl,k,Z=Z,fullW0=W0)
      fit = variance_adjust(fit, bin=bin)
      my_individual_analysis(fit, Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, k=k, do_rle=TRUE, do_projection=TRUE)
    }
  }

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

  if(do_ruvinv)
  {
    fit = RUVinv(Y,X,ctl,Z=Z)
    fit = variance_adjust(fit, bin=bin)
    # No RLE plots with RUVinv
    my_individual_analysis(fit, Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, do_projection=TRUE)
  }

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

  if(do_ruvrinv)
  {
    fit = RUVrinv(Y,X,ctl,Z=Z)
    fit = variance_adjust(fit, bin=bin)
    my_individual_analysis(fit, Y, pctl, samplecoloring, genecoloring, genetexts, genesymbols, topcount_thresholds, geneinfo, topN, rankset=rankset, do_projection=TRUE)
  }

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

  if(do_pptable)
  {
    fullfit = RUVinv(Y,X,ctl,Z=Z,inputcheck=FALSE)
    my_projplottable(fullfit, c(0,kset), factorset, genecoloring, genetexts, genesymbols)
  }
}

######################################################################################
## General Info Routine
######################################################################################

my_general_info = function(m, n, nc, Y, Yc, X, Z, ctl, W_fullctl, W_ctl, samplecoloring, genecoloring, sampletexts, genetexts, samplesymbols, genesymbols)
{
  # Display m, n, n_c
  cat("m = ", m, ", n = ", n, ", n_c = ", nc, "\n")
  # Make scree plots, adjusted for Z
  my_ruv_scree(Y, Z=Z)
  my_ruv_scree(Yc, Z=Z)
  # Make RLE plots.  Don't adjust for Z, since we'll make Z-adjusted RLE plots in the "Unadjusted" section anyway.
  my_ruv_rle(Y, XZW=NULL, coloring=samplecoloring)
  my_ruv_rle(Yc, XZW=NULL, coloring=samplecoloring)
  # Make canonical correlation plots.
  my_ruv_cancorplot(X, W_fullctl, W2=W_ctl, col=c(1, genecoloring[ctl][1]))
  # Make PC plots.
  my_ruv_pcgridplot(W_fullctl, W2=W_ctl, kset=c(1,2,3), coloring=samplecoloring, symbols=samplesymbols, texts=sampletexts)
  # Make alpha plots.  Same routine as above, but now we just pass in alpha values.
  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)
}

######################################################################################
## Individual Analysis Routine
######################################################################################

my_individual_analysis <-
function(fit, 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)
{
  if (do_rle)
  {
      my_ruv_rle(Y, XZW=cbind(fit$X,fit$Z,fit$W), coloring=samplecoloring)	
      my_ruv_rle(Y[,fit$ctl,drop=FALSE], XZW=cbind(fit$X,fit$Z,fit$W), coloring=samplecoloring)
  }
  if (do_pval)
  {
      my_ruv_pvalhist(fit$p)
      my_ruv_pvalhist(fit$p.ebayes)
      my_ruv_pvalhist(fit$p.rsvar)
      my_ruv_pvalhist(fit$p.rsvar.ebayes)
      my_ruv_pvalhist(fit$p.evar)
      my_ruv_pvalplot(fit$p, coloring=genecoloring, symbols=genesymbols, texts=genetexts)
      my_ruv_pvalplot(fit$p.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts)
      my_ruv_pvalplot(fit$p.rsvar, coloring=genecoloring, symbols=genesymbols, texts=genetexts)
      my_ruv_pvalplot(fit$p.rsvar.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts)
      my_ruv_pvalplot(fit$p.evar, coloring=genecoloring, symbols=genesymbols, texts=genetexts)
      my_ruv_pvalplot(fit$p, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)
      my_ruv_pvalplot(fit$p.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)
      my_ruv_pvalplot(fit$p.rsvar, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)
      my_ruv_pvalplot(fit$p.rsvar.ebayes, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)
      my_ruv_pvalplot(fit$p.evar, coloring=genecoloring, symbols=genesymbols, texts=genetexts, logged=TRUE)
  }
  if (do_ctlpval)
  {
      my_ruv_pvalhist(fit$p[fit$ctl])
      my_ruv_pvalhist(fit$p.ebayes[fit$ctl])
      my_ruv_pvalhist(fit$p.rsvar[fit$ctl])
      my_ruv_pvalhist(fit$p.rsvar.ebayes[fit$ctl])
      my_ruv_pvalhist(fit$p.evar[fit$ctl])
      my_ruv_pvalplot(fit$p[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])
      my_ruv_pvalplot(fit$p.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])    
      my_ruv_pvalplot(fit$p.rsvar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])
      my_ruv_pvalplot(fit$p.rsvar.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])
      my_ruv_pvalplot(fit$p.evar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl])
      my_ruv_pvalplot(fit$p[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)
      my_ruv_pvalplot(fit$p.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)
      my_ruv_pvalplot(fit$p.rsvar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)
      my_ruv_pvalplot(fit$p.rsvar.ebayes[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)
      my_ruv_pvalplot(fit$p.evar[fit$ctl], coloring=genecoloring[fit$ctl], symbols=genesymbols[fit$ctl], texts=genetexts[fit$ctl], logged=TRUE)
  }
  if (do_projection)
  {
      my_ruv_projectionplot(fit$byx, fit, fit$bwx, coloring=genecoloring, symbols=genesymbols, texts=genetexts)
      my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p, rankset=rankset)
      my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.ebayes, rankset=rankset)
      my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.rsvar, rankset=rankset)
      my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.rsvar.ebayes, rankset=rankset)
      my_ruv_topprojectionplot(fit$byx, fit, fit$bwx, fit$p.evar, rankset=rankset)
}
  if (do_variance)
  {
      my_ruv_varianceplot(fit, fit$varbetahat, coloring=genecoloring, symbols=genesymbols, texts=genetexts, trans=my_fourthroot)
      my_ruv_topvarianceplot(fit, fit$p, fit$varbetahat, trans=my_fourthroot, rankset=rankset)
      my_ruv_topvarianceplot(fit, fit$p.ebayes, fit$varbetahat.ebayes, trans=my_fourthroot, rankset=rankset)
      my_ruv_topvarianceplot(fit, fit$p.rsvar, fit$varbetahat.rsvar, trans=my_fourthroot, rankset=rankset)
      my_ruv_topvarianceplot(fit, fit$p.rsvar.ebayes, fit$varbetahat.rsvar.ebayes, trans=my_fourthroot, rankset=rankset)
      my_ruv_topvarianceplot(fit, fit$p.evar, fit$varbetahat.evar, trans=my_fourthroot, rankset=rankset)
  }
  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)
      cat("Number of Top-ranked Positive Controls\n")
      cat("standard\n")
      print(rankcount)
      cat("ebayes\n")
      print(rankcount.ebayes)
      cat("rsvar\n")
      print(rankcount.rsvar)
      cat("rsvar ebayes\n")
      print(rankcount.rsvar.ebayes)
      cat("evar\n")
      print(rankcount.evar)
    }
  }
  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)
    cat("Top-ranked Genes\n")
    cat("standard\n")
    print(ranktable)
    cat("ebayes\n")
    print(ranktable.ebayes)
    cat("rsvar\n")
    print(ranktable.rsvar)
    cat("rsvar ebayes\n")
    print(ranktable.rsvar.ebayes)
    cat("evar\n")
    print(ranktable.evar)
  }
}

######################################################################################
## Combined Analysis Routine
######################################################################################

my_combined_analysis <-
function(counttables)
{
  if (!is.null(counttables))
  {
    my_ruv_toprankcountplot(counttables$countmat, counttables$kset, counttables$thresholds)
    my_ruv_toprankcountplot(counttables$countmat.ebayes, counttables$kset, counttables$thresholds)
    my_ruv_toprankcountplot(counttables$countmat.rsvar, counttables$kset, counttables$thresholds)
    my_ruv_toprankcountplot(counttables$countmat.rsvar.ebayes, counttables$kset, counttables$thresholds)
    my_ruv_toprankcountplot(counttables$countmat.evar, counttables$kset, counttables$thresholds)
  }
}

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

my_projplottable <- function(fullfit, kset, factorset, genecoloring, genetexts, genesymbols)
{
  for (i in 1:length(kset)) for (j in 1:length(factorset))
  {
    my_ruv_factor_projectionplot(fullfit, kset[i], factorset[j], coloring=genecoloring, symbols=genesymbols, texts=genetexts)
  }
}