nsc <-
function (x, y, n.threshold = 30, offset.percent = 50, prior = NULL, remove.zeros = TRUE) 
{
  this.call <- match.call()
  Y <- model.matrix(~factor(y) - 1, data = list(y = y))
  
  xtest <- x
  ytest <- y
  
  
  n.class <- table(y)
  if (min(n.class) == 1) {
    stop(warning("Warning: a class contains only 1 sample"))
  }

  n <- sum(n.class)
  ntest <- ncol(xtest)
  K <- length(prior)
  p <- nrow(x)
  
  dimnames(Y) <- list(NULL, names(n.class))
  centroids <- scale(x %*% Y, FALSE, n.class)  ## WMEAN.G
  
  xdif <- x - centroids %*% t(Y)
  sd <- (xdif^2) %*% rep(1/(n - K), n)
  sd <- drop(sqrt(sd))  #WSD.POOLED
  offset <- quantile(sd, offset.percent/100)
  sd <- sd + offset

  centroid.overall <- drop(x %*% rep(1/n, n))   ## WMEAN
  
  se.scale <- sqrt(1/n.class - 1/n)  # mk
  
  delta <- (centroids - centroid.overall)/sd
  delta <- scale(delta, FALSE, se.scale)    ##dik
  
  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/(se.scale))
    delta.shrunk <- t(t(delta.shrunk) * as.numeric(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 <- safe.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, prob = prob[, , seq(n.threshold)], 
                 centroids = centroids, centroid.overall = centroid.overall, 
                 sd = sd, threshold = threshold[seq(n.threshold)], nonzero = nonzero[seq(n.threshold)], 
                 se.scale = se.scale, call = this.call, prior = prior, offset = offset)
  if (!is.null(ytest)) 
    object$errors <- errors
  #class(object) <- "nsc"
  object
}