--- a
+++ b/predict.voomNSC.R
@@ -0,0 +1,89 @@
+predict.voomNSC =
+function (fit, newdata, threshold = NULL, prior = NULL) 
+{
+  if (is.null(prior)) prior = fit$prior
+  if (is.null(threshold)) threshold = fit$opt.threshold
+  vm = voomGSD(data.train = fit$counts, data.test = newdata, group = fit$conditions, norm.method = fit$normalization)
+  x = vm$TestExp
+  
+  
+# soft.shrink(...)
+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
+  }
+
+# diag.disc(...)
+diag.disc =
+  function (x, centroids, prior, weight) 
+  {
+    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)
+  }
+
+# safe.exp(...)
+safe.exp = function (x) 
+{
+  xx = sign(x) * pmin(abs(x), 500)
+  return(exp(xx))
+}
+
+#soft.max(...)
+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
+  }
+
+  sd = fit$weightedSD.pooled
+  centroid.overall = fit$weightedMean
+  centroids = fit$weightedMean.C
+  se.scale = fit$se.scale
+  delta = fit$delta
+ 
+  delta.shrunk = soft.shrink(delta, threshold)
+  delta.shrunk = t(t(delta.shrunk) * as.numeric(se.scale))  
+  posid = drop(abs(delta.shrunk) %*% rep(1, length(prior))) > 0
+  dd = diag.disc((x - fit$weightedMean)/(fit$weightedSD.pooled), delta.shrunk, prior, posid)
+  softmax(dd)
+}
\ No newline at end of file