source('http://hivdb.stanford.edu/pages/published_analysis/phenotypeCompJAID2005/statScripts/ROC_lib/functions.R')
source(Rurl('read_data.R'))
source(Rurl('cutoffs.R'))

protein <- 'RT'

ROCs <- list()

pdf('ROC.pdf')

for (which.subset in c(1,2)) {
  for (drug in drugs) {
    .eval(paste('ROCs$', drug, '.', patternnames[[which.subset]], ' <- list()', sep=''))
    for (assay in assays){
      cur.table <- get.data(assay, protein, type='wild')
      wild <- log(cur.table$Fold[(cur.table$Drug == drug)])
      cur.table <- get.data(assay, protein, type='mut')
      mut <- log(cur.table$Fold[(cur.table$Drug == drug) & (cur.table$subset == which.subset)])
      cur.ROC.table <- list()

      cur.ROC.table$wild.ecdf <- ecdf(wild)
      cur.ROC.table$mut.ecdf <- ecdf(mut)

      cur.ROC.table$ROC.emp <- get.ROC.emp(wild, mut)
      cur.ROC.table$AUC.emp <- integrate(approxfun(cur.ROC.table$ROC.emp$x, cur.ROC.table$ROC.emp$y), 0, 1)$value

      cur.ROC.table$ROC.dens <- get.ROC.dens(wild, mut)
      cur.ROC.table$AUC.dens <- integrate(approxfun(cur.ROC.table$ROC.dens$x, cur.ROC.table$ROC.dens$y), 0, 1)$value

      .eval(paste('ROCs$', drug, '.', patternnames[[which.subset]], '$', assay, ' <- cur.ROC.table', sep=''))
    }

    upper.table <- get.ROC.table(drug, which.subset)

    virco.wild.ecdf <- upper.table$virco$wild.ecdf
    virologic.wild.ecdf <- upper.table$virologic$wild.ecdf

    virco.mut.ecdf <- upper.table$virco$mut.ecdf
    virologic.mut.ecdf <- upper.table$virologic$mut.ecdf

    virco.cutoffs <- .eval(paste('cutoffs$', drug, '$virco', sep=''))
    virologic.cutoffs <- .eval(paste('cutoffs$', drug, '$virologic', sep=''))

    test <- ((drug != 'D4T') || (which.subset != 2))
    if (test) {
      plot(upper.table$virco$ROC.emp, type='n', col='blue', lwd=2, xlab='False Positive Rate', ylab='True Positive Rate', main=paste(drug, ':', group[[which.subset]]))

      xval <- seq(-4,4,0.01)
      lines(1 - virologic.wild.ecdf(xval), 1 - virologic.mut.ecdf(xval), type='l', lwd=2, col='red')
      lines(1 - virco.wild.ecdf(xval), 1 - virco.mut.ecdf(xval), type='l', lwd=2, col='blue')

      symbol <- c(21:23)
      for (i in 1:3) {
        points(1 - virco.wild.ecdf(virco.cutoffs[i]), 1 - virco.mut.ecdf(virco.cutoffs[i]), pch=symbol[i], bg='yellow', cex=1.5)
        points(1 - virologic.wild.ecdf(virologic.cutoffs[i]), 1 - virologic.mut.ecdf(virologic.cutoffs[i]), pch=symbol[i], bg='yellow', cex=1.5)
      }
    }

  }
}

dev.off()