.packageName <- "pamr"
 balanced.folds <- function(y, nfolds = min(min(table(y)), 10)) {
   totals <- table(y)
   fmax <- max(totals)
   nfolds <- min(nfolds, fmax)     
                                        # makes no sense to have more folds than the max class size
   folds <- as.list(seq(nfolds))
   yids <- split(seq(y), y)        
                                        # nice we to get the ids in a list, split by class
###Make a big matrix, with enough rows to get in all the folds per class
   bigmat <- matrix(NA, ceiling(fmax/nfolds) * nfolds, length(totals))
   for(i in seq(totals)) {
     bigmat[seq(totals[i]), i] <- sample(yids[[i]])
   }
   smallmat <- matrix(bigmat, nrow = nfolds)       # reshape the matrix
### Now do a clever sort to mix up the NAs
   smallmat <- permute.rows(t(smallmat))   ### Now a clever unlisting
                                        # the "clever" unlist doesn't work when there are no NAs
                                        #       apply(smallmat, 2, function(x)
                                        #        x[!is.na(x)])
   res <-vector("list", nfolds)
   for(j in 1:nfolds) {
     jj <- !is.na(smallmat[, j])
     res[[j]] <- smallmat[jj, j]
   }
   return(res)
 }
descendants <- function(m,k){
  ## the done object indicates what rows of m were used
  done <- k
  if (m[k,1] < 0)
    left <- -m[k,1]
  else {
    junk <- descendants(m, m[k,1])
    left <- junk[[1]]
    done <- c(done, junk[[2]])
  }
  if (m[k,2] < 0)
    right <- -m[k,2]
  else {
    junk <- descendants(m, m[k,2])
    right <- junk[[1]]
    done <- c(done, junk[[2]])
  } 
  return(list(c(left, right), done))
}
diag.disc <-function(x, centroids, prior, weight) {
### Computes the class discriminant functions assuming scaled x and centroids
  if(! missing(weight)) {
    posid <- (weight > 0)
    if(any(posid)) {
      weight <- sqrt(weight[posid])
      centroids <- centroids[posid,  , drop = FALSE] * weight
      x <- x[posid,  , drop = FALSE] * weight
    }
    else {
      mat <- outer(rep(1, ncol(x)), log(prior), "*")
      dimnames(mat) <- list(NULL, dimnames(centroids)[[2]])
      return(mat)
    }
  }
  dd <- t(x) %*% centroids
  dd0 <- drop(rep(1, nrow(centroids)) %*% (centroids^2))/2 - log(prior)
  names(dd0) <- NULL
  scale(dd, dd0, FALSE)
}
enlist <-function(...) {
  result <- list(...)
  if((nargs() == 1) & is.character(n <- result[[1]])) {
    result <- as.list(seq(n))
    names(result) <- n
    for(i in n)
      result[[i]] <- get(i)
  }
  else {
    junk <- sys.call()
    n <- NULL
    for(i in junk[-1])
      n <- c(n, deparse(i))
    if(!is.null(n2 <- names(result))) {
      which <- n2 != ""
      n[which] <- n2[which]
    }
    names(result) <- n
  }
  result
}
error.nsc <-function(object) {
###Computes the roc curve for a nsc model
  yhat <- object$yhat
  y <- object$y
  ny <- table(y)
  errors <- matrix(0, length(object$threshold), length(ny))
  Y <- data.matrix(yhat) != unclass(y)
  yind <- model.matrix( ~ factor(y) - 1, data = list(y = y))
  errors <- t(t(yind) %*% Y)
  apply(errors, 2, mean)
}
nsc <-function(x, y, xtest = x, ytest = NULL, threshold = NULL, n.threshold = 30, 
        hetero=NULL,
        scale.sd = TRUE, threshold.scale = NULL, se.scale = NULL, offset.percent=50, prior = table(y)/length(y), remove.zeros = TRUE, sign.contrast="both")
{
        this.call <- match.call()

        n.class <- table(y)
if(min(n.class)==1){stop("Error: each class must have >1 sample")}

        norm.cent <-NULL
        if(!is.null(hetero)){
           norm.cent <-apply(x[,y==hetero],1,mean)
           x <-abs(t(scale(t(x),center=norm.cent,scale=FALSE)))
          if(!missing(xtest)){xtest <-abs(t(scale(t(xtest),center=norm.cent,scale=FALSE)))}
        }

        n <- sum(n.class)
        ntest <- ncol(xtest)
        K <- length(prior)
        p <- nrow(x)
        if(missing(xtest))
                ytest <- y
        Y <- model.matrix( ~ factor(y) - 1, data = list(y = y))
        dimnames(Y) <- list(NULL, names(n.class))
        centroids <- scale(x %*% Y, FALSE, n.class)
        sd <- rep(1, p)
        if(scale.sd) {
                xdif <- x - centroids %*% t(Y)
                sd <- (xdif^2) %*% rep(1/(n - K), n)
                sd <- drop(sqrt(sd))
                offset  <- quantile(sd, offset.percent/100)
                sd <- sd + offset
        }
        centroid.overall <- drop(x %*% rep(1/n, n))
        if(is.null(threshold.scale)) {
                threshold.scale <- rep(1, K)
                names(threshold.scale) <- names(n.class)
        }
### Now make an adjustment for the sample sizes in the "t" ratios
        if(is.null(se.scale))
                se.scale <- sqrt(1/n.class - 1/n)
        delta <- (centroids - centroid.overall)/sd
        delta <- scale(delta, FALSE, threshold.scale * se.scale)    

 if(sign.contrast=="positive"){delta <- delta*(delta>0)}
  if(sign.contrast=="negative"){delta <- delta*(delta<0)}



        #allows differential shrinkage
        if(!is.null(threshold)) {
                n.threshold <- length(threshold)
        }
        else {
                threshold <- seq(0, max(abs(delta)), length = n.threshold)
        }
        nonzero <- seq(n.threshold)
        errors <- threshold
        yhat <- as.list(seq(n.threshold))
        prob <- array(0, c(ntest, K, n.threshold))
        for(ii in 1:n.threshold) {
                cat(ii)
                delta.shrunk <- soft.shrink(delta, threshold[ii])
                delta.shrunk <- scale(delta.shrunk, FALSE, 1/(threshold.scale * 
                        se.scale))
                nonzero[ii] <- attr(delta.shrunk, "nonzero")
                posid <- drop(abs(delta.shrunk) %*% rep(1, K)) > 0
                dd <- diag.disc((xtest - centroid.overall)/sd, delta.shrunk, 
                        prior, weight = posid)
                yhat[[ii]] <- softmax(dd)
                dd <- exp(dd)
                prob[,  , ii] <- dd/drop(dd %*% rep(1, K))
                if(!is.null(ytest)) {
                        errors[ii] <- sum(yhat[[ii]] != ytest)
                }
        }
        thresh.names <- format(round(threshold, 3))
        names(yhat) <- thresh.names
        attr(yhat, "row.names") <- paste(seq(ntest))
        class(yhat) <- "data.frame"
        if(remove.zeros)
                n.threshold <- match(0, nonzero, n.threshold)
        dimnames(prob) <- list(paste(seq(ntest)), names(n.class), thresh.names)
        object <- list(y = ytest, yhat = yhat[, seq(n.threshold)], prob = 
                       prob[,  , seq(n.threshold)], centroids=centroids, centroid.overall=centroid.overall, sd=sd, 
                       threshold = threshold[seq(n.threshold)], nonzero = nonzero[seq(
                n.threshold)], threshold.scale=threshold.scale, se.scale=se.scale, call = this.call, hetero=hetero,
                       norm.cent=norm.cent,
                prior=prior, offset=offset, sign.contrast=sign.contrast)
        if(!is.null(ytest))
                object$errors <- errors[seq(n.threshold)]
        class(object) <- "nsc"
        object
}
nsccv <-
function(x, y, nfold = min(table(y)), folds = balanced.folds(y), threshold =
        NULL, threshold.scale = NULL, prior, object, ...)
{
        this.call <- match.call()
        n <- length(y)
        if(is.null(folds)) {
                folds <- split(sample(1:n), rep(1:nfold, length = n))
        }
        else nfold <- length(folds)
        if(missing(prior)) {
                if(missing(object))
                        prior <- table(y)/n
                else prior <- object$prior
        }
        if(missing(threshold)) {
                if(missing(object))
                        stop("Must either supply threshold argument, or an nsc object"
                                )
                else {
                        threshold <- object$threshold
                        threshold.scale <- object$threshold.scale
                        se.scale <- object$se.scale
                }
        }
        n.threshold <- length(threshold)        ### Set up the data structures
        yhat <- rep(list(y), n.threshold)
        names(yhat) <- paste(seq(n.threshold))
        yhat <- data.frame(yhat)
        n.class <- table(y)
        prob <- array(1, c(n, length(n.class), n.threshold))
        size <- double(n.threshold)
        hetero <-object$hetero
        for(ii in 1:nfold) {
                cat("Fold", ii, ":")
                a <- nsc(x[,  - folds[[ii]]], y[ - folds[[ii]]], x[, folds[[ii
                        ]], drop = FALSE], threshold = threshold, threshold.scale
                         = threshold.scale, se.scale = se.scale, prior = prior,
                          hetero=hetero,
                        ..., remove.zeros = FALSE)
                size <- size + a$nonzero
                prob[folds[[ii]],  ,  ] <- a$prob
                yhat[folds[[ii]],  ] <- a$yhat
                cat("\n")
        }
        if(missing(object))
                size <- round(size/nfold)
        else size <- object$nonzero
        error <- rep(NA, n.threshold)
        loglik <- error
        for(i in 1:n.threshold) {
                error[i] <- sum(yhat[, i] != y)/n
                loglik[i] <- sum(log(prob[,  , i][cbind(seq(1, n), unclass(y))]))/                        n
        }
obj<- list(threshold=threshold, error=error, loglik=loglik,size=size, yhat=yhat,y=y,prob=prob,folds=folds,
                call = this.call)
        class(obj) <- "nsccv"
        obj
}

pamr.adaptthresh <- function(object, ntries = 10, reduction.factor = 0.9, full.out = FALSE) {
  errors <- error.nsc(object)
  threshold <- object$threshold   
### Remove all but the first leading zero errors
  ifirst <- match(TRUE, object$errors > 0, FALSE)
  if (!ifirst)
    stop("Zero training error throughout!")
  else {
    ifirst <- max(ifirst, 1)
    threshold <- threshold[seq(ifirst, length(threshold))]
  }
### initialization
  tscales <- object$threshold.scale
  all.errors <- matrix(0, ntries + 1, length(tscales),
                       dimnames = list(NULL, names(tscales)))
  all.scales <- all.errors
  all.objects <- as.list(seq(ntries + 1))
  rocs <- double(ntries + 1)
  all.scales[1,  ] <- tscales
  all.errors[1,  ] <- errors
  rocs[1] <- roc.nsc(object)      # integrated size^(1/4)*error
  cat("Initial errors:", format(round(errors, 5)), "Roc",
      format(round(rocs[1], 5)), "\n")
  for (i in seq(ntries)) {
    cat("Update", i, "\n")
    j <- rev(order(errors))[1]      # identify the largest error
    tscales[j] <- tscales[j] * reduction.factor     
                                        # and reduce its scale
    all.scales[i + 1,  ] <- tscales/min(tscales)    # and renormalize
    iobject <- update(object, threshold = threshold, 
                      threshold.scale = all.scales[i + 1,  ], remove.zeros = 
                      FALSE)
    all.errors[i + 1,  ] <- errors <- error.nsc(iobject)
    rocs[i + 1] <- roc.nsc(iobject)
    cat("\nErrors", format(round(errors, 5)), "Roc",
        format(round(rocs[i + 1], 5)), "\n")
  }
  j <- order(rocs)[1]     # identify the scales with the smallest "roc"
  opt.scale <- all.scales[j,  ]
  if (full.out)
    list(errors = all.errors, scales = all.scales, rocs = rocs, 
         opt.scale = opt.scale)
  else
    opt.scale
}

pamr.batchadjust <- function(data) {
  if (is.null(data$batchlabels)) {
    stop("batch labels are not in data object")
  }
  lab <- data$batchlabels
  dd <- model.matrix( ~ factor(lab) - 1)
  data$x <- data$x - misreg.simple(dd, data$x)
  data
}


misreg.simple <- function(Y, x) {
###Y is a indicator response matrix
  nax <- is.na(x)
  nsamples <- (!nax)%*%Y
  x[nax] <- 0
  xsum <- x%*%Y
  xbar <- xsum/nsamples
  xbar %*% t(Y)
}
pamr.confusion <- function(fit, threshold, extra = TRUE) {
  ii <- (1:length(fit$threshold))[fit$threshold > threshold]
  ii <- ii[1]
  predicted <- fit$yhat[, ii]
  if (is.null(fit$newy)) {
    true <- fit$y[fit$sample.subset]
  }
  else {
    true <- fit$newy[fit$sample.subset]
  }
  tt <- table(true, predicted)
  if (extra) {
    tt1 <- tt
    diag(tt1) <- 0
    tt <- cbind(tt, apply(tt1, 1, sum)/apply(tt, 1, sum))
    dimnames(tt)[[2]][ncol(tt)] <- "Class Error rate"
    print(tt)
    cat(c("Overall error rate=", round(sum(tt1)/sum(tt), 3)),
        fill= TRUE)
  }
  if (!extra) {
    return(tt)
  }
}
## Cube root transformation for Affy chips
pamr.cube.root  <- function(x) {
  return(sign(x) * abs(x)^{1/3})
}

pamr.cv <-
function(fit, data, nfold = min(table(data$y)), folds = balanced.folds(data$y),...)
{
        x <- data$x[fit$gene.subset, fit$sample.subset]
        if(is.null(fit$newy)) {
                y <- factor(data$y[fit$sample.subset])
        }
        else {
                y <- factor(data$newy[fit$sample.subset])
        }
        this.call <- match.call()
        junk <- nsccv(x, y, object = fit, ...)
        junk$call <- this.call
        junk$newy <- fit$newy
        junk$sample.subset <- fit$sample.subset
        return(junk)
}

pamr.from.excel <- function(file, ncols, sample.labels = FALSE, batch.labels = FALSE) {
  d <- scan(file, sep = "\t", what = "")
  dd <- matrix(d, ncol = ncols, byrow = TRUE)
  samplelabels <- NULL
  batchlabels <- NULL
  ii <- 1
  if(sample.labels) {
    samplelabels <- dd[1,  - (1:2)]
    ii <- ii + 1
  }
  if(batch.labels & !sample.labels) {
    batchlabels <- dd[1,  - (1:2)]
    ii <- ii + 1
  }
  if(batch.labels & sample.labels) {
    batchlabels <- dd[2,  - (1:2)]
    ii <- ii + 1
  }
  y <- dd[ii,  - (1:2)]
  geneid <- dd[ - (1:ii), 1]
  genenames <- dd[ - (1:ii), 2]
  x <- matrix(as.numeric(as.character(dd[ - (1:ii),  - (1:2)])), ncol = 
              ncols - 2)
  cat("",fill=TRUE)
  cat(c("Read in ", nrow(x), "genes"),fill=TRUE)
  cat(c("Read in ", ncol(x), "samples"),fill=TRUE)
  if(sample.labels){cat(c("Read in ", length(samplelabels), "sample labels"),fill=TRUE)}
  if(batch.labels){cat(c("Read in ", length(batchlabels), "batch labels"),fill=TRUE)}
  cat("",fill=TRUE)
  cat("Make sure these figures are correct!!", fill=TRUE)
  cat("",fill=TRUE)
  
  return(list(x = x, y = y, genenames = genenames, geneid = geneid, 
              samplelabels = samplelabels, batchlabels = batchlabels))
}

pamr.geneplot <- function(fit, data, threshold) {
  par(pch = 1, col = 1)
  geneid <- data$geneid
  if(is.null(geneid)) {
    geneid <- as.character(1:nrow(data$x))
  }
  if(is.null(fit$newy)) {
    y <- factor(data$y[fit$sample.subset])
  }
  else {
    y <- factor(fit$newy[fit$sample.subset])
  }
  x <- data$x[fit$gene.subset, fit$sample.subset]
  geneid <- geneid[fit$gene.subset]
  nc <- length(unique(y))
  aa <- pamr.predict(fit, x, threshold = threshold, type = "nonzero")
  cen <- pamr.predict(fit, x, threshold = threshold, type = "cen")
  d <- (cen - fit$centroid.overall)[aa,  ]/fit$sd[aa]
  oo <- order( - apply(abs(d), 1, max))
  aa <- aa[oo]
  ngenes <- length(aa)
  o <- order(y)
  xx <- x[aa, o]
  geneid <- geneid[aa]
  nc <- length(unique(y))
  nn <- c(0, cumsum(table(y)))
  nrow <- trunc(sqrt(ngenes)) + 1
  ncol <- trunc(sqrt(ngenes)) + 1
  if(nrow * (ncol - 1) >= ngenes) {
    ncol <- ncol - 1
  }
  par(mfrow = c(nrow, ncol))
  for(i in 1:ngenes) {
    plot(1:ncol(xx), xx[i,  ], type = "n", xlab = "sample", ylab = 
         "expression", axes = FALSE)
    box()
    axis(2)
    for(j in 1:nc) {
      j1 <- nn[j] + 1
      j2 <- nn[j] + table(y)[j]
      points(j1:j2, xx[i, j1:j2], col = j + 1)
    }
    title(main = as.character(geneid[i]))
    for(j in 1:(nc - 1)) {
      abline(v = cumsum(table(y))[j] + 0.5, lty = 2)
    }
    if(i == 1) {
      h <- c(0, table(y))
      for(j in 2:(nc + 1)) {
        text(sum(h[1:(j - 1)]) + 0.5 * h[j], max(xx[i,  
                                                    ]), label = levels(y)[j - 1], col = j)
      }
    }
  }
  par(mfrow = c(1, 1))
}


pamr.knnimpute <- function(data, k = 10) {
  x <- data$x
  N <- dim(x)
  p <- N[2]
  
  N <- N[1]
        col.nas  <- apply(x, 2, is.na)
  if ((sum(col.nas) == N) > 0) {
    stop("Error: A column has all missing values!")
  }
  nas <- is.na(drop(x %*% rep(1, p)))
  xcomplete <- x[!nas,  ]
  xbad <- x[nas,,drop=FALSE ]
  xnas <- is.na(xbad)
  xbadhat <- xbad
  cat(nrow(xbad), fill = TRUE)
  for(i in seq(nrow(xbad))) {
    cat(i, fill = TRUE)
    xinas <- xnas[i,  ]
    xbadhat[i,  ] <- nnmiss(xcomplete, xbad[i,  ], xinas, K = k)
  }
  x[nas,  ] <- xbadhat
  data2 <-data
  data2$x <-x
  return(data2)
}

nnmiss <- function(x, xmiss, ismiss, K = 1) {
  xd <- scale(x, xmiss, FALSE)[, !ismiss]
  dd <- drop(xd^2 %*% rep(1, ncol(xd)))
  od <- order(dd)[seq(K)]
  xmiss[ismiss] <- drop(rep(1/K, K) %*% x[od, ismiss, drop = FALSE])
  xmiss
}

pamr.listgenes <- function (fit, data, threshold, genenames = FALSE)  {
  if (is.null(fit$newy)) {
    y <- factor(data$y[fit$sample.subset])
  }
  if (!is.null(fit$newy)) {
    y <- factor(fit$newy[fit$sample.subset])
  }
  x <- data$x[fit$gene.subset, fit$sample.subset]
  if (genenames) {
    gnames <- data$genenames[fit$gene.subset]
  }
  if (!genenames) {
    gnames <- NULL
  }
  geneid <- data$geneid[fit$gene.subset]
  nc <- length(unique(y))
  aa <- pamr.predict(fit, x, threshold = threshold, type = "nonzero")
  cen <- pamr.predict(fit, x, threshold = threshold, type = "cen")
  d <- (cen - fit$centroid.overall)[aa, ]/fit$sd[aa]
  
  oo <- order(-apply(abs(d), 1, max))
  d <- round(d, 4)
  g <- gnames[aa]
  g1 <- geneid[aa]
  if (is.null(gnames)) {
    gnhdr <- NULL
  }
  if (!is.null(gnames)) {
    gnhdr <- "name"
  }
  options(width = 500)
  schdr <- paste(dimnames(table(y))$y, "score", sep = " ")
  res <- cbind(as.character(g1), g, d)[oo, ]
  dimnames(res) <- list(NULL, c("id", gnhdr, schdr))
  print(res, quote = FALSE)
}

pamr.makeclasses <- function(data,  sort.by.class = FALSE, ...) {
  require(cluster)
  as.matrix.dist <- function (x)  {
    size <- attr(x, "Size")
    df <- matrix(0, size, size)
    df[row(df) > col(df)] <- x
    df <- df + t(df)
    labels <- attr(x, "Labels")
    dimnames(df) <- if (is.null(labels)) 
      list(1:size, 1:size)
    else list(labels, labels)
    df
  }
  as.dist <- function (m, diag = FALSE, upper = FALSE) {
    m <- as.matrix(m)
    retval <- m[row(m) > col(m)]
    attributes(retval) <- NULL
    if (!is.null(rownames(m))) 
      attr(retval, "Labels") <- rownames(m)
    else if (!is.null(colnames(m))) 
      attr(retval, "Labels") <- colnames(m)
    attr(retval, "Size") <- nrow(m)
    attr(retval, "Diag") <- diag
    attr(retval, "Upper") <- upper
    attr(retval, "call") <- match.call()
    class(retval) <- "dist"
    retval
  }
  
  if(!is.null(data$samplelabels)) {
    labs <- data$samplelabels
  }
  if(!is.null(data$samplelabels) & !is.null(data$y)) {
    labs <- paste(data$y, labs)
  }
  if(is.null(data$samplelabels)) {
    labs <- 1:ncol(data$x)
  }
  par(col = 1, cex = 1)
  d <- dist(t(data$x))
  dd <- as.matrix.dist(d)
  if(sort.by.class) {
    tt <- table(data$y)
    nc <- length(tt)
    for(i in 1:nc) {
      o <- data$y == names(tt[i])
      d1 <- max(dd[o, o])
      d2 <- min(dd[o, !o])
      fac <- ((0.2 + (0.7 * i)/nc) * d2)/d1
      dd[o, o] <- dd[o, o] * fac
    }
  }
  hc <- hclust(as.dist(dd), ...)
  plot(hc, labels = labs)
  aa <- vector("list", 100)
  go <- TRUE
  i <- 0
  while(go & i < 100) {
    go <- FALSE
    i <- i + 1
    print(c("Identify class", i))
    par(pch = as.character(i), col = 4)
    aa[[i]] <- locator(type = "p")
    if(!is.null(aa[[i]])) {
      go <- TRUE
    }
  }
  nclus <- i - 1
  res <- vector("list", nclus)
  for(i in 1:nclus) {
    res[i] <- aa[i]
  }
  hdelta <- 1
  clus <- vector("list", nclus)
  for(j in 1:nclus) {
    for(jj in 1:length(res[[j]]$x)) {
      r <- c(res[[j]]$x[jj], res[[j]]$y[jj])
      d <- abs(hc$hei - r[2])
      o <- rank(d)
      ncomp <- 5
      oo <- (1:length(o))[o < ncomp + 1 & d < hdelta]
      if(length(oo) == 0) {
        stop(
             "1 Ambigious selection; try pamr.makeclasses again"
             )
      }
      ncomp2 <- length(oo)
      good <- rep(FALSE, ncomp2)
      ordpos <- match(1:length(hc$ord), hc$ord)
      nodes <- vector("list", ncomp2)
      for(ii in 1:ncomp2) {
        ooo <- descendants(hc$mer, oo[ii])[[2]]
        o4 <- as.vector(hc$mer[ooo,  ])
        nodes[[ii]] <- -1 * o4[o4 < 0]
        op <- ordpos[nodes[[ii]]]
        if(r[1] > min(op) & r[1] < max(op)) {
          good[ii] <- TRUE
        }
      }
                                        #browser()
      if(sum(good) != 1) {
        stop(
             "2 Ambigious selection; try pamr.makeclasses again"
             )
      }
                                        #browser()
      ii2 <- (1:ncomp2)[good]
      clus[[j]] <- c(clus[[j]], nodes[[ii2]])
    }
  }
  newy <- rep(NA, ncol(data$x))
  temp <- NULL
  for(i in 1:nclus) {
    clus[[i]] <- unique(clus[[i]])
  }
  for(i in 1:nclus) {
    temp <- c(temp, clus[[i]])
  }
  if(length(unique(temp)) < length(temp)) {
    stop("Clusters overlap; try pamr.makeclasses again")
  }
  for(i in 1:nclus) {
    newy[clus[[i]]] <- i
  }
  labs2 <- as.character(newy)
  labs2[labs2 == "NA"] <- ""
  par(col = 1, cex = 1)
  plot(hc, labels = labs2)
  return(as.factor(newy))
}

pamr.menu <- function(data) {
  done <- FALSE
  junk.train <- NULL
  junk.results <- NULL
  while(!done) {
    cat("", fill = TRUE)
    switch(menu(c("pamr.train", "pamr.cv", "pamr.plotcv", 
                  "pamr.plotcen", "pamr.confusion", 
                  "pamr.plotcvprob", "pamr.geneplot", 
                  "pamr.listgenes", 
                  "pamr.train with heterogeneity analysis", 
                  "Exit")),
           junk.train <- pamr.train(data),
           {
             if(is.null(junk.train)) {
               cat("Error: need to run pamr.train first", 
                   fill = TRUE)
             }
             if(!is.null(junk.train)) {
               junk.results <- pamr.cv(junk.train, data)
             }
           }
           ,
           {
             if(is.null(junk.results)) {
               cat("Error: need to run pamr.cv first", fill
                   = TRUE)
             }
             if(!is.null(junk.results)) {
               pamr.plotcv(junk.results)
             }
           }
           ,
           {
             if(is.null(junk.train)) {
               cat("Error: need to run pamr.train first", 
                   fill = TRUE)
             }
             if(!is.null(junk.train)) {
               cat("threshold?")
               threshold <- scan("", nlines = 1)
               pamr.plotcen(junk.train, data, threshold = 
                            threshold)
             }
           }
           ,
           {
             if(is.null(junk.results)) {
               cat("Error: need to run pamr.cv first", fill
                   = TRUE)
             }
             if(!is.null(junk.results)) {
               cat("threshold?")
               threshold <- scan("", nlines = 1)
               pamr.confusion(junk.results, threshold = 
                              threshold)
             }
           }
           ,
           {
             if(is.null(junk.results)) {
               cat("Error: need to run pamr.cv first", fill
                   = TRUE)
             }
             if(!is.null(junk.results)) {
               cat("threshold?")
               threshold <- scan("", nlines = 1)
               pamr.plotcvprob(junk.results, data, threshold
                               = threshold)
             }
           }
           ,
           {
             if(is.null(junk.train)) {
               cat("Error: need to run pamr.train first", 
                   fill = TRUE)
             }
             if(!is.null(junk.train)) {
               cat("threshold?")
               threshold <- scan("", nlines = 1)
               pamr.geneplot(junk.train, data, threshold = 
                             threshold)
             }
           }
           ,
           {
             if(is.null(junk.train)) {
               cat("Error: need to run pamr.train first", 
                   fill = TRUE)
             }
             if(!is.null(junk.train)) {
               cat("threshold?")
               threshold <- scan("", nlines = 1)
               pamr.listgenes(junk.train, data, threshold = 
                              threshold)
             }
           }
           ,
           {
             junkk.train <- NULL
             cat("Normal class?", fill = TRUE)
             normal <- scan("", nlines = 1, what = "")
             junk.train <- pamr.train(data, hetero = normal)
           }
           ,
           done <- TRUE)
  }
  cat("Done\n")
}

pamr.pairscore <-function(x, pair.ind=NULL) {
}

pamr.plotcen <- function(fit, data, threshold) {
  genenames <- data$genenames[fit$gene.subset]
  x <- data$x[fit$gene.subset, fit$sample.subset]
  clabs <- colnames(fit$centroids)
  scen <- pamr.predict(fit, data$x, threshold = threshold, type = "cent")
  dif <- scen - fit$centroid.overall
  nc <- length(unique(fit$y))
  o <- drop(abs(dif) %*% rep(1, nc)) > 0
  d <- dif[o,  ]
  nd <- sum(o)
  genenames <- genenames[o]
  xx <- x[o,  ]
  oo <- order(apply(abs(d), 1, max))
  d <- d[oo,  ]
  genenames <- genenames[oo]
  par(mar = c(1, 5, 1, 1), col = 1)
  plot(rep(2, nd) + d[, 1], 1:nd, xlim = c(0, 2*nc+1), ylim = c(1, nd + 3), 
       type = "n", xlab = "", ylab = "", axes = FALSE)
  box()
  abline(h = seq(nd), lty = 3, col = 7)
  jj <- rep(0, nd)
  for(j in 1:nc) {
    segments(jj + 2 * j, seq(nd), jj + 2 * j + d[, j], seq(nd), col
             = j + 1, lwd = 4)
    lines(c(2 * j, 2 * j), c(1, nd), col = j + 1)
    text(2 * j, nd + 2, label = clabs[j], col = j + 1)
  }
  g <- substring(genenames, 1, 20)
  text(rep(0, nd), seq(nd), label = g, cex = 0.4, adj = 0, col = 1)
}

pamr.plotcen <- function(fit, data, threshold) {
  genenames <- data$genenames[fit$gene.subset]
  x <- data$x[fit$gene.subset, fit$sample.subset]
  clabs <- colnames(fit$centroids)
  scen <- pamr.predict(fit, data$x, threshold = threshold, type = "cent")
  dif <- (scen - fit$centroid.overall)/fit$sd
  nc <- length(unique(fit$y))
  o <- drop(abs(dif) %*% rep(1, nc)) > 0
  d <- dif[o,  ]
  nd <- sum(o)
  genenames <- genenames[o]
  xx <- x[o,  ]
  oo <- order(apply(abs(d), 1, max))
  d <- d[oo,  ]
  genenames <- genenames[oo]
  par(mar = c(1, 5, 1, 1), col = 1)
  plot(rep(2, nd) + d[, 1], 1:nd, xlim = c(0, 2*nc+1), ylim = c(1, nd + 3), 
       type = "n", xlab = "", ylab = "", axes = FALSE)
  box()
  abline(h = seq(nd), lty = 3, col = 7)
  jj <- rep(0, nd)
  for(j in 1:nc) {
    segments(jj + 2 * j, seq(nd), jj + 2 * j + d[, j], seq(nd), col
             = j + 1, lwd = 4)
    lines(c(2 * j, 2 * j), c(1, nd), col = j + 1)
    text(2 * j, nd + 2, label = clabs[j], col = j + 1)
  }
  g <- substring(genenames, 1, 20)
  text(rep(0, nd), seq(nd), label = g, cex = 0.4, adj = 0, col = 1)
}
pamr.plotcv <- function(fit) {
  par(mar = c(5, 5, 5, 1))
  par(mfrow = c(2, 1))
  n <- nrow(fit$yhat)
  y <- fit$y
  if(!is.null(fit$newy)) {
    y <- fit$newy[fit$sample.subset]
  }
  nc <- length(table(y))
  nfolds <- length(fit$folds)
  err <- matrix(NA, ncol = ncol(fit$yhat), nrow = nfolds)
  temp <- matrix(y, ncol = ncol(fit$yhat), nrow = n)
  ni <- rep(NA, nfolds)
  for(i in 1:nfolds) {
    ii <- fit$folds[[i]]
    ni[i] <- length(fit$folds[[i]])
    err[i,  ] <- apply(temp[ii,  ] != fit$yhat[ii,  ], 2, sum)/ni[i]
  }
  se <- sqrt(apply(err, 2, var)/nfolds)
  plot(fit$threshold, fit$error, ylim = c(-0.1, 0.8), xlab = 
       "Value of threshold  ", ylab = "Misclassification Error", type
       = "n", yaxt = "n")
  axis(3, at = fit$threshold, lab = paste(fit$size), srt = 90, adj = 0)
  mtext("Number of genes", 3, 4, cex = 1.2)
  axis(2, at = c(0, 0.2, 0.4, 0.6, 0.8))
  lines(fit$threshold, fit$error, col = 2)
  o <- fit$err == min(fit$err)
  points(fit$threshold[o], fit$error[o], pch = "x")
  error.bars(fit$threshold, fit$err - se, fit$err + se)
  err2 <- matrix(NA, nrow = length(unique(y)), ncol = length(fit$threshold
                                                 ))
  for(i in 1:(length(fit$threshold) - 1)) {
    s <- pamr.confusion(fit, fit$threshold[i], extra = FALSE)
    diag(s) <- 0
    err2[, i] <- apply(s, 1, sum)/table(y)
  }
  plot(fit$threshold, err2[1,  ], ylim = c(-0.1, 1.1), xlab = 
       "Value of threshold ", ylab = "Misclassification Error", type
       = "n", yaxt = "n")
  axis(3, at = fit$threshold, lab = paste(fit$size), srt = 90, adj = 0)     
                                        #       mtext("Number of genes", 3, 4,cex=1.2)
  axis(2, at = c(0, 0.2, 0.4, 0.6, 0.8))
  for(i in 1:nrow(err2)) {
    lines(fit$threshold, err2[i,  ], col = i + 1)
  }
  legend(0, 0.9, dimnames(table(y))[[1]], col = (2:(nc + 1)), lty = 1)
  par(mfrow = c(1, 1))
}

error.bars <-function(x, upper, lower, width = 0.02, ...) {
  xlim <- range(x)
  barw <- diff(xlim) * width
  segments(x, upper, x, lower, ...)
  segments(x - barw, upper, x + barw, upper, ...)
  segments(x - barw, lower, x + barw, lower, ...)
  range(upper, lower)
}

pamr.plotcvprob <- function(fit, data, threshold) {
  par(pch = 1)
  ii <- (1:length(fit$threshold))[fit$threshold > threshold]
  ii <- ii[1]
  ss <- data$samplelabels
  pp <- fit$prob[,  , ii]
  if(is.null(fit$newy)) {
    y <- fit$y[fit$sample.subset]
  }
  if(!is.null(fit$newy)) {
    y <- fit$newy[fit$sample.subset]
  }
  o <- order(y)
  y <- y[o]
  if(!is.null(ss)) {
    ss <- ss[o]
  }
  ppp <- pp[o,  ]
  n <- nrow(ppp)
  nc <- length(unique(y))
  par(cex = 1)
  plot(1:n, ppp[, 2], type = "n", xlab = "sample", ylab = 
       "cross-validated probabilities", ylim = c(0, 1.2), axes = FALSE)
  axis(1)
  axis(2, labels = c("0.0", "0.2", "0.4", "0.6", "0.8", "1.0", ""))
  axis(4)
  for(j in 1:nc) {
    points(1:n, ppp[, j], col = j + 1)
  }
  for(j in 1:(nc - 1)) {
    abline(v = cumsum(table(y))[j] + 0.5, lty = 2)
  }
  h <- c(0, table(y))
  for(j in 2:(nc + 1)) {
    text(sum(h[1:(j - 1)]) + 0.5 * h[j], 1.02, label = levels(y)[j - 
                                                 1], col = j)
  }
  abline(h = 1)
  if(!is.null(ss)) {
    text(1:length(ss), 1.1, labels = ss, srt = 90, cex = 0.7)
  }
  ##if(!is.null(ss)){axis(3,labels=ss,at=1:length(ss),srt=90)}
}

pamr.predict <-  function(fit, newx, threshold, type = c("class", "posterior", "centroid", "nonzero"), 
                          prior = fit$prior,  threshold.scale = fit$
                          threshold.scale) {
  norm.cen <- fit$norm.cen
  if(!is.null(norm.cen)) {
    newx <- abs(t(scale(t(newx), center = norm.cen, scale = FALSE)))
  }
  type <- match.arg(type)
  sd <- fit$sd
  centroid.overall <- fit$centroid.overall
  centroids <- fit$centroids
  se.scale <- fit$se.scale
  delta <- scale((centroids - centroid.overall)/sd, FALSE, threshold.scale * 
                 se.scale)

  if(fit$sign.contrast=="positive"){delta <- delta*(delta>0)}
  if(fit$sign.contrast=="negative"){delta <- delta*(delta<0)}


  delta.shrunk <- scale(soft.shrink(delta, threshold), FALSE, 1/(
                                                                 threshold.scale * se.scale))
  posid <- drop(abs(delta.shrunk) %*% rep(1, length(prior))) > 0
  if(!match(type, c("centroid", "nonzero"), FALSE))
    dd <- diag.disc((newx - centroid.overall)/sd, delta.shrunk, 
                    prior, posid)
  switch(type,
         class = softmax(dd),
         posterior = {
           dd <- exp(dd)
           dd/drop(dd %*% rep(1, length(prior)))
         }
         ,
         centroid = centroid.overall + delta.shrunk * sd,
         nonzero = {
           nz <- drop(abs(delta.shrunk) %*% rep(1, ncol(centroids)
                                                )) > 0
           seq(nz)[nz]
         }
         )
}

pamr.predictmany <- function(fit, newx, threshold=fit$threshold,
                             prior = fit$prior,  threshold.scale = fit$threshold.scale,
                             ...) {
  prob <-array(NA,c(length(prior),ncol(newx),length(threshold)))
  predclass <-matrix(NA,nrow=ncol(newx),ncol=length(threshold))
  
  for(i in 1:length(threshold)){
    prob[,,i] <-pamr.predict(fit,newx,threshold=threshold[i],type="posterior",...)
    predclass[,i] <-pamr.predict(fit,newx,threshold=threshold[i],type="class",...)
  }
  
  predclass <-matrix(levels(fit$y)[predclass],ncol=length((threshold)))

  return(list(prob=prob,predclass=predclass))
}













pamr.to.excel <- function(data, file, trace = TRUE) {
  if(is.null(data$x) | is.null(data$y) | is.null(data$genenames) | 
     is.null(data$geneid)) {
    stop("Invalid format for input data")
  }
  n <- nrow(data$x)
  p <- ncol(data$x)
  row1 <- paste("", "", sep = "\t")
  if(!is.null(data$samplelabels)) {
    for(j in 1:p) {
      row1 <- paste(row1, data$samplelabels[j], sep = "\t")
    }
    write(row1, file = file, append = FALSE)
  }
  row2 <- paste("", "", sep = "\t")
  if(!is.null(data$batchlabels)) {
    for(j in 1:p) {
      row2 <- paste(row2, data$batchlabels[j], sep = "\t")
    }
    write(row2, file = file, append = TRUE)
  }
  row3 <- paste("", "", sep = "\t")
  for(j in 1:p) {
    row3 <- paste(row3, data$y[j], sep = "\t")
  }
  write(row3, file = file, append = TRUE)
  for(i in 1:n) {
    if(trace) {
      cat(c("writing row number", i), fill = TRUE)
    }
    xx <- paste(data$genenames[i], data$geneid[i], sep = "\t")
    for(j in 1:ncol(data$x)) {
      xx <- paste(xx, data$x[i, j], sep = "\t")
    }
    write(xx, file = file, append = TRUE)
  }
  return()
}
pamr.train <-

function(data, gene.subset=1:nrow(data$x), sample.subset=1:ncol(data$x),
         threshold = NULL, n.threshold = 30,
        scale.sd = TRUE, threshold.scale = NULL, se.scale = NULL, offset.percent = 50, hetero=NULL,
         prior = NULL,  remove.zeros = TRUE, sign.contrast="both")

{
        this.call <- match.call()

        if(is.null(prior))
          {prior <- table(data$y[sample.subset])/length(data$y[sample.subset])
           prior <- prior[prior!=0]
        }

        if(!is.null(sample.subset) & !is.null(data$newy)) {
           stop("Can't have both newy present in data object, and sample.subset specified"
                        )
                     }
       if(is.null(sample.subset)){sample.subset <-1:ncol(data$x)}

        if(is.null(data$y) & is.null(data$newy)) {
                stop("must have either y or newy present in data object")
        }
        if(is.null(data$newy)) {
                y <- data$y
        }
        if(!is.null(data$newy)) {
                y <- data$newy
                sample.subset <- (1:ncol(data$x))[!is.na(y)]
                print("Using classes `newy' from data object")
        }
        junk <- nsc(data$x[gene.subset, sample.subset], factor(y[sample.subset]), 
          offset.percent=offset.percent,  threshold = threshold, hetero=hetero,
          n.threshold = n.threshold,  scale.sd= scale.sd, threshold.scale=threshold.scale,
           se.scale= se.scale, prior=prior, remove.zeros=remove.zeros,
            sign.contrast=sign.contrast)

        junk$call <- this.call
        junk$gene.subset <- gene.subset
        junk$sample.subset <- sample.subset
        junk$newy <- data$newy
        return(junk)
}
pamr.xl.compute.offset <- function(data, offset.percent=50, prior=prior){
  x <- data$x
  y <- data$y
  n.class <- table(y)
  if(min(n.class)==1){stop("Error: each class must have >1 sample")}
  norm.cent <-NULL
  n <- sum(n.class)
  xtest <- x
  ntest <- ncol(xtest)
  K <- length(prior)
  p <- nrow(x)
  Y <- model.matrix( ~ factor(y) - 1, data = list(y = y))
  dimnames(Y) <- list(NULL, names(n.class))
  centroids <- scale(x %*% Y, FALSE, n.class)
  sd <- rep(1, p)
  xdif <- x - centroids %*% t(Y)
  sd <- (xdif^2) %*% rep(1/(n - K), n)
  sd <- drop(sqrt(sd))
  offset  <- quantile(sd, offset.percent/100)
  return(offset)
}

pamr.xl.get.offset  <- function() {
  if (exists("x.train")) {
    return (x.train$offset)
  } else {
    return (pamr.xl.compute.offset(pamr.xl.data, offset.percent=pamr.xl.training.parameters$offset.percent,
                                   prior=pamr.xl.training.parameters$prior))
  }
}

pamr.xl.derive.adjusted.prior  <- function(prior, data) {
  s  <- pamr.xl.get.sample.prior(data)
  temp <- prior - s
  if (sum(temp*temp) < pamr.xl.training.parameters$epsilon) {
    return (list (prior=s, prior.name="Sample Prior"))
  } else {
    s  <-  pamr.xl.get.uniform.prior(data)
    temp  <- prior - s
    if (sum(temp*temp) < pamr.xl.training.parameters$epsilon) {
      return (list (prior=s, prior.name="Uniform Prior"))
    } else {
      return (list (prior=prior, prior.name="Custom Prior"))      
    }
  }
}

pamr.xl.get.default.training.parameters <- function(data) {
  return (list(offset.percent=50, prior=pamr.xl.get.sample.prior(data), prior.name="Sample Prior", sign.contrast="both", epsilon=1e-7))
}

## Return the uniform prior on class labels
pamr.xl.get.uniform.prior  <- function(data) {
  w <- table(data$y)
  n  <- length(w)
  return(rep(1.0/n, n))
}

## Return the sample proportion prior on class labels
pamr.xl.get.sample.prior  <- function(data) {
  w <- table(data$y)
  return(w/sum(w))
}


pamr.xl.process.data <- function() {
  res <- list(x=pamr.xl.raw.data, y=pamr.xl.class.labels, genenames=pamr.xl.gene.names, 
              geneid=pamr.xl.gene.ids, samplelabels=pamr.xl.sample.labels,
              batchlabels=pamr.xl.batch.labels)
  
  if (pamr.xl.data.has.missing.values) {
    res <- pamr.knnimpute(res, k = pamr.xl.knn.neighbors)
  }
  return(res)
}

pamr.xl.compute.cv.confusion  <- function (fit, cv.results, threshold) {
  threshold.rank  <- which(rank(abs(cv.results$threshold - threshold))==1)
  t.threshold  <- cv.results$threshold[threshold.rank]
  true  <- fit$y
  predicted  <- cv.results$yhat[, threshold.rank]
  tt <- table(true, predicted)
  tt1 <- tt
   diag(tt1) <- 0
  tt <- cbind(tt, apply(tt1, 1, sum)/apply(tt, 1, sum))
  dimnames(tt)[[2]][ncol(tt)] <- "Class Error rate"
  overall.err  <- round(sum(tt1)/sum(tt), 3)
  return(list(confusion.matrix=tt, overall.error=overall.err, threshold=round(t.threshold, 5)))
 }
pamr.xl.compute.confusion  <- function (fit, threshold) {
  ii <- (1:length(fit$threshold))[fit$threshold > threshold]
  ii <- ii[1]
  predicted <- fit$yhat[, ii]
  if(is.null(fit$newy)) {
    true <- fit$y[fit$sample.subset]
  }
  else {
    true <- fit$newy[fit$sample.subset]
  }
  tt <- table(true, predicted)
  tt1 <- tt
  diag(tt1) <- 0
  tt <- cbind(tt, apply(tt1, 1, sum)/apply(tt, 1, sum))
  dimnames(tt)[[2]][ncol(tt)] <- "Class Error rate"
  overall.err  <- round(sum(tt1)/sum(tt), 3)
  return(list(confusion.matrix=tt, overall.error=overall.err))
}
pamr.xl.is.a.subset  <- function(x, y) {
  if (nlevels(factor(x)) == nlevels(factor(c(x, y[!is.na(y)])))) {
    return (1)  # True
  } else {
    return (0)  # False
  }
}
pamr.xl.listgenes.compute  <- function (fit, data, threshold, genenames = FALSE) {
  if (is.null(fit$newy)) {
    y <- factor(data$y[fit$sample.subset])
  }
  if (!is.null(fit$newy)) {
    y <- factor(fit$newy[fit$sample.subset])
  }
  x <- data$x[fit$gene.subset, fit$sample.subset]
  if (genenames) {
    gnames <- data$genenames[fit$gene.subset]
  }
  if (!genenames) {
    gnames <- NULL
  }
  geneid <- data$geneid[fit$gene.subset]
  nc <- length(unique(y))
  aa <- pamr.predict(fit, x, threshold = threshold, type = "nonzero")
  cen <- pamr.predict(fit, x, threshold = threshold, type = "cen")
  d <- (cen - fit$centroid.overall)[aa, ]/fit$sd[aa]
  oo <- order(-apply(abs(d), 1, max))
  d <- round(d, 4)
  g <- gnames[aa]
  g1 <- geneid[aa]
  if (is.null(gnames)) {
    gnhdr <- NULL
  }
  if (!is.null(gnames)) {
    gnhdr <- "name"
  }
  options(width = 500)
  schdr <- paste(dimnames(table(y))$y, "score", sep = " ")
  res <- cbind(as.character(g1), g, d)[oo, ]
  dimnames(res) <- list(NULL, c("id", gnhdr, schdr))
  return(list(gene.headings = dimnames(res)[[2]],
              gene.ids = res[ , 2],   # This was switched with gene.names. 
              gene.names = res[ , 1],
              gene.scores = res[ , -(1:2)]))
  ##print(res, quote = FALSE)
}
pamr.xl.plot.test.probs.compute  <- function(fit, new.x, newx.classes, missing.class.label, 
	threshold, sample.labels=NULL) {
  predicted.probs  <- pamr.xl.predict.test.probs(x.train, new.x, threshold=threshold)
  training.classes  <- levels(factor(fit$y))
  py  <- pamr.xl.predict.test.class.only(x.train, new.x, threshold=threshold)
  order.classes  <- order(newx.classes)
  pp  <- predicted.probs[, order.classes]
  actual.classes <- newx.classes[order.classes]
  actual.classes[is.na(actual.classes)] <- missing.class.label
  ny  <- py$predicted[order.classes]
  n  <- length(ny)
  ss  <- sample.labels
  if (!is.null(ss)) {
    ss  <- ss[order.classes]
  }
  
  return (list(x = 1:n,
               y = t(pp),
               x.label = "Sample",
               y.label = "Predicted Test Probabilities",
               y.names = levels(factor(fit$y)),
               y.lines = cumsum(table(actual.classes)) + 0.5,
               x.dummy = vector(length=2, mode="numeric"),
               y.dummy = vector(length=2, mode="numeric"),
               panel.names = levels(factor(actual.classes)),
               x.names = ss))
}  



pamr.xl.plot.training.error.compute  <- function(trained.object) {
  return (list(x = trained.object$threshold,
               y = trained.object$errors/length(trained.object$y),
               y.ytop = trained.object$nonzero,
               x.label = "Threshold",
               y.label = "Training Error"))
}
pamr.xl.plotcen.compute  <- function(fit, data, threshold) {
  genenames <- data$genenames[fit$gene.subset]
  x <- data$x[fit$gene.subset, fit$sample.subset]
  clabs <- colnames(fit$centroids)
  scen <- pamr.predict(fit, data$x, threshold = threshold, type = "cent")
  dif <- scen - fit$centroid.overall
  nc <- length(unique(fit$y))
  o <- drop(abs(dif) %*% rep(1, nc)) > 0
  d <- dif[o,  ]
  nd <- sum(o)
  genenames <- genenames[o]
  xx <- x[o,  ]
  oo <- order(apply(abs(d), 1, max))
  d <- d[oo,  ]
  genenames <- genenames[oo]
  win.metafile()
  plot.title=paste("Centroid Plot( Threshold =", threshold, ")")
  par(mar = c(1, 5, 1, 1), col = 1)
  plot(rep(2, nd) + d[, 1], 1:nd, xlim = c(0, 2*nc+1), ylim = c(1, nd + 3), 
       type = "n", xlab = "", ylab = "", axes = FALSE, main=plot.title)
  box()
  abline(h = seq(nd), lty = 3, col = 7)
  jj <- rep(0, nd)
  for(j in 1:nc) {
    segments(jj + 2 * j, seq(nd), jj + 2 * j + d[, j], seq(nd), col
             = j + 1, lwd = 4)
    lines(c(2 * j, 2 * j), c(1, nd), col = j + 1)
    text(2 * j, nd + 2, label = clabs[j], col = j + 1)
  }
  g <- substring(genenames, 1, 20)
  text(rep(0, nd), seq(nd), label = g, cex = 0.4, adj = 0, col = 1)
  dev.off()
#  pamr.plot.y <<- matrix(d, nrow=dim(d)[1])
#  pamr.plot.x <<- seq(nd)
#  pamr.plot.seriesnames <<- dimnames(d)[[2]]
#  pamr.plot.genenames <<- genenames

  return(TRUE)
}
pamr.xl.plotcv.compute  <- function(aa) {
  n <- nrow(aa$yhat)
  y <- aa$y
  if(!is.null(aa$newy)) {
    y <- aa$newy[aa$sample.subset]
  }
  nc <- length(table(y))
  nfolds <- length(aa$folds)
  err <- matrix(NA, ncol = ncol(aa$yhat), nrow = nfolds)
  temp <- matrix(y, ncol = ncol(aa$yhat), nrow = n)
  ni <- rep(NA, nfolds)
  for(i in 1:nfolds) {
    ii <- aa$folds[[i]]
    ni[i] <- length(aa$folds[[i]])
    err[i,  ] <- apply(temp[ii,  ] != aa$yhat[ii,  ], 2, sum)/ni[i]
  }
  se <- sqrt(apply(err, 2, var)/nfolds)

  err2 <- matrix(NA, nrow = length(unique(y)), ncol = length(aa$threshold)-1)
  for(i in 1:(length(aa$threshold) - 1)) {
    s <- pamr.confusion(aa, aa$threshold[i], extra = FALSE)
    diag(s) <- 0
    err2[, i] <- apply(s, 1, sum)/table(y)
  }

  return (list(x = aa$threshold,
               y = aa$error,
               x.label = "Threshold",
               y.label = "Misclasiffication Error",
               y.se = se,
               y.ytop = aa$size,
               cv.err = t(err2),
               cv.legend = dimnames(table(y))[[1]]))
               
}
pamr.xl.plotcvprob.compute  <- function(aa, data, threshold) {
  ii <- (1:length(aa$threshold))[aa$threshold > threshold]
  ii <- ii[1]
  ss <- data$samplelabels
  pp <- aa$prob[,  , ii]
  if(is.null(aa$newy)) {
    y <- aa$y[aa$sample.subset]
  }
  if(!is.null(aa$newy)) {
    y <- aa$newy[aa$sample.subset]
  }
  o <- order(y)
  y <- y[o]
  if(!is.null(ss)) {
    ss <- ss[o]
  }
  ppp <- pp[o,  ]
  n <- nrow(ppp)
  nc <- length(unique(y))


#  axis(2, labels = c("0.0", "0.2", "0.4", "0.6", "0.8", "1.0", ""))
#  if (!is.null(ss)) {
#    pamr.plot.x.names <<- ss
#  }

  return (list(x = 1:n,
               y = ppp,
               x.label = "Sample",
               y.label = "CV Probabilities",
               y.names = levels(y),
               y.lines = cumsum(table(data$y)),
               x.dummy = vector(length=2, mode="numeric"),
               y.dummy = vector(length=2, mode="numeric"),
               x.names = ss))
  
#   for(j in 1:nc) {
#     points(1:n, ppp[, j], col = j + 1)
#   }
#   for(j in 1:(nc - 1)) {
#     abline(v = cumsum(table(y))[j] + 0.5, lty = 2)
#   }
#   h <- c(0, table(y))
#   for(j in 2:(nc + 1)) {
#     text(sum(h[1:(j - 1)]) + 0.5 * h[j], 1.02, label = levels(y)[j - 
#                                                  1], col = j)
#   }
#   abline(h = 1)
#   if(!is.null(ss)) {
#     text(1:length(ss), 1.1, labels = ss, srt = 90, cex = 0.7)
#   }
  ##if(!is.null(ss)){axis(3,labels=ss,at=1:length(ss),srt=90)}
}
pamr.xl.predict.test.class<- function(fit, newx, threshold, test.class.labels) {
  predicted  <- pamr.predict(fit, newx, threshold, type="class")
  return(list(confusion.matrix=table(test.class.labels, predicted), predicted=as.vector(predicted)))
}

pamr.xl.predict.test.class.only  <- function(fit, newx, threshold) {
  return(list(predicted=as.vector(pamr.predict(fit, newx, threshold, type="class"))))
}
pamr.xl.predict.test.class.only  <- function(fit, newx, threshold) {
  return(list(predicted=as.vector(pamr.predict(fit, newx, threshold, type="class"))))
}
pamr.xl.predict.test.probs  <- function(fit, newx, threshold) {
  predicted  <- pamr.predict(fit, newx, threshold, type="posterior")
  return(t(predicted))
}

pamr.xl.test.data.impute  <- function(x, k) {
  N <- dim(x)
  p <- N[2]
  N <- N[1]
  col.nas  <- apply(x, 2, is.na)
  if ((sum(col.nas) == N) > 0) {
    stop("Error: A column has all missing values!")
  }
  
  nas <- is.na(drop(x %*% rep(1, p)))
  xcomplete <- x[!nas,  ]
  xbad <- x[nas,,drop=FALSE ]
  xnas <- is.na(xbad)
  xbadhat <- xbad
  cat(nrow(xbad), fill = TRUE)
  for(i in seq(nrow(xbad))) {
    cat(i, fill = TRUE)
    xinas <- xnas[i,  ]
    xbadhat[i,  ] <- nnmiss(xcomplete, xbad[i,  ], xinas, K = k)
  }
  x[nas,  ] <- xbadhat
  return(x)
}
pamr.xl.test.errors.compute  <- function(fit, newx, newx.classes, threshold=fit$threshold,
                                         prior = fit$prior,  threshold.scale = fit$threshold.scale,
                                         ...) {
  n  <- length(which(!is.na(newx.classes)))
## Note: n is assumed to be nonzero! Check before calling!
  actual.classes  <- newx.classes
  prediction.errs  <- vector(mode="numeric", length=length(threshold))
  
  for(i in 1:length(threshold)){
    t <- pamr.predict(fit,newx,threshold=threshold[i],type="class",...)
    prediction.errs[i]  <- length(which(t != actual.classes)) / n
  }
  
  return(list(x=threshold, y=prediction.errs, x.label="Threshold", y.label="Test Error", ))
  
}
pamr.xl.transform.class.labels  <- function(x) {
  y  <- x
  y[is.na(y)]  <- " "
  return(y)
}

pamr.xl.transform.data <- function(data) {

  if (pamr.xl.take.cube.root) {
    data$x = pamr.cube.root(data$x)
  }

  if (pamr.xl.batch.labels.present) {
    data <- pamr.batchadjust(data)
  }

  if (pamr.xl.center.columns && pamr.xl.scale.columns) {
    data$x = scale(data$x, center=TRUE, scale=TRUE)
  } else if (pamr.xl.center.columns) {
    data$x = scale(data$x, center=TRUE, scale=FALSE)
  } else if (pamr.xl.scale.columns) {
    data$x = scale(data$x, center=FALSE, scale=TRUE)
  }

  return (data)
}

pamr.xl.transform.test.data <- function(test.x) {
  res <- test.x
  if (pamr.xl.take.cube.root) {
    res = pamr.cube.root(res)
  }

  if (pamr.xl.center.columns && pamr.xl.scale.columns) {
    res = scale(res, center=TRUE, scale=TRUE)
  } else if (pamr.xl.center.columns) {
    res = scale(res, center=TRUE, scale=FALSE)
  } else if (pamr.xl.scale.columns) {
    res = scale(res, center=FALSE, scale=TRUE)
  }

  return (res)
}

permute.rows <-function(x)
{
        dd <- dim(x)
        n <- dd[1]
        p <- dd[2]
        mm <- runif(length(x)) + rep(seq(n) * 10, rep(p, n))
        matrix(t(x)[order(mm)], n, p, byrow = TRUE)
}

print.nsc <- function(x, ...) {
  cat("Call:\n")
  dput(x$call)
  mat <- rbind(threshold = format(round(x$threshold, 3)), nonzero = 
               format(trunc(x$nonzero)), errors = x$errors)
  dimnames(mat) <- list(dimnames(mat)[[1]], paste(1:ncol(mat)))
  print(t(mat), quote = FALSE)
  invisible()
}
 print.nsccv <-function(x, ...) {
   cat("Call:\n")
   dput(x$call)
   mat <- rbind(threshold = format(round(x$threshold, 3)), nonzero = 
                format(trunc(x$size)), errors = trunc(x$error * nrow(
                                              x$yhat)))
   dimnames(mat) <- list(dimnames(mat)[[1]], paste(1:ncol(mat)))
   print(t(mat), quote = FALSE)
   invisible()
 }
roc.nsc <-function(object) {
###Computes the roc curve for a nsc model
  nonzero <- object$nonzero^(1/4)
  errors <- object$errors
  if(is.null(errors))
    stop("No errors component")
  n <- length(errors)
  heights <- (errors[1:(n - 1)] + errors[2:n])/2
  bases <- diff(nonzero)
  area <- sum((nonzero[-1] + nonzero[-n]) * heights * bases) /
    (-2 * diff(range(nonzero)))
  area
}
softmax <-function(x, gap = FALSE) {
  d <- dim(x)
  maxdist <- x[, 1]
  pclass <- rep(1, d[1])
  for(i in seq(2, d[2])) {
    l <- x[, i] > maxdist
    pclass[l] <- i
    maxdist[l] <- x[l, i]
  }
  dd <- dimnames(x)[[2]]
  if(gap) {
    x <- abs(maxdist - x)
    x[cbind(seq(d[1]), pclass)] <- drop(x %*% rep(1, d[2]))
    gaps <- do.call("pmin", data.frame(x))
  }
  pclass <- if(is.null(dd) || !length(dd))
    pclass
  else
    factor(pclass, levels = seq(d[2]), labels = dd)
  if(gap)
    list(class = pclass, gaps = gaps)
  else
    pclass
}
soft.shrink <-function(delta, threshold) {
  dif <- abs(delta) - threshold
  delta <- sign(delta) * dif * (dif > 0)
  nonzero <- sum(drop((dif > 0) %*% rep(1, ncol(delta))) > 0)
  attr(delta, "nonzero") <- nonzero
  delta
}
