fisherRRSVD <- function(X,E,lambda=1,scale.by.sv=FALSE) {

  if (missing(X)) {
    cat("Usage: V <- fisherRRSVD(X,E,lambda)
      X is nSamples x nComponents
      T is nSamples x 1 column vector of integer class labels
      lambda is regularization amount.  0 means no regularization
      V is matrix of discriminate directions in columns

  from \"A Flexible and Efficient Algorithm for Regularized Fisher
  Discriminant Analysis\", by Z. Zhang, G. Dai, and M. I. Jordan, 2009.
  In Proceedings of the  Machine Learning and Knowledge Discovery in
  Databases: European Conference (ECML PKDD), Bled, Slovenia,
  pages 632-647, 2009.
  www.cs.berkeley.edu/~jordan/papers/zhang-dai-jordan-ecml09.pdf\n")
    return(invisible(NULL))
  }

  n <- nrow(X)
  p <- ncol(X)

  if (is.null(dim(E)) || (length(dim(E)) == 2 && ncol(E)==1))
    E <- makeIndicatorVars(E)
  PiDiag <- colSums(E)

  X <- X - matrix(colMeans(X),n,p,byrow=TRUE)
  E <- E / matrix(sqrt(PiDiag),n,ncol(E),byrow=TRUE)
  XE <- t(X) %*% E
  
  if (n < p) {
    ##use (14)
    invP <- solve(X %*% t(X) + lambda * diag(1,n), E)
    G <- t(X) %*% invP
    R <- t(XE) %*% G
    ##R 		= (XE'*X') * invP; 			#this might be faster; needs profiling;
  } else {
    ##use (13)
    G <- solve(t(X) %*% X + lambda * diag(1,p), XE)
    R <- t(XE) %*% G
  }

  eigenResult <- eigen(0.5 * (R+t(R)))
  V <- eigenResult$vectors
  gam <- eigenResult$values
  lastOne <- ncol(V)
  A <- G %*% V[,-lastOne,drop=FALSE]
  if (scale.by.sv) {
    scalby <- sqrt(gam[-lastOne])
    scalby <- scalby
    A <- A / matrix(scalby,nrow(A),ncol(A),byrow=TRUE)
  }
  A / matrix(sqrt(colSums(A^2)),nrow(A),ncol(A),byrow=TRUE)
}

######################################################################
######################################################################
######################################################################
######################################################################

makeIndicatorVars <- function(Y) {
  if (!is.matrix(Y))
    Y <- matrix(Y)
  classes <- unique(Y)
  N <- nrow(Y)
  K <- length(classes)
  logicalMatrix <- (matrix(Y,N,K) == matrix(classes,N,K,byrow=TRUE))
  mode(logicalMatrix) <- "numeric" ## to convert to numbers 0, 1
  logicalMatrix
}