#
#                 Generalized Linear Discriminant Analysis
#
#              Copyright(c) Haifeng Li <hli@cs.ucr.edu> 2004
#                           All rights reserved
#

library(class)

orth <- function( u )
{
  v <- u
  v[,1] <- v[,1]/sqrt(t(v[,1])%*%v[,1])
  if ( ncol(u) > 1 )
  {
    for ( i in 2 : ncol(u) )
    {
      for ( j in 1 : (i-1) )
      {
        v[,i] <- v[,i] - (t(u[,i])%*%v[,j])*v[,j]
      }
      v[,i] <- v[,i]/sqrt(t(v[,i])%*%v[,i])
    }
  }
  v
}

glda <- function( x, g )
{
  n <- nrow(x)
  p <- ncol(x)
  c <- nlevels(g)
  prior <- table(g)/n
  clusters <- as.numeric(levels(g))

  group.means <- tapply(x, list(rep(g, p), col(x)), mean)
  total.mean <- apply(group.means, 2, weighted.mean, prior)

  M <- matrix(0, nrow=c, ncol=p)
  for ( i in 1:c )
  {
    M[i,] <- sqrt(prior[i])*(group.means[i,]-total.mean)
  }
  M <- t(M)

  X <- matrix(0, nrow=n, ncol=p)
  for ( i in 1:n )
  {
    X[i,] <- (x[i,]-total.mean)/sqrt(n)
  }
  X <- t(X)

  XSVD <- svd(X)

  dinv <- matrix(0, length(XSVD$d), length(XSVD$d))
  for ( i in 1 : length(XSVD$d) )
    if ( abs(XSVD$d[i]) > 1E-5 ) dinv[i,i] <- 1/XSVD$d[i]

  A <- t(XSVD$u) %*% M
  A <- dinv %*% A
  A <- XSVD$u %*% A
  ASVD <- svd(A)
  W <- t(XSVD$u) %*% ASVD$u[,1:(c-1)]
  W <- dinv %*% W
  W <- XSVD$u %*% W
  W <- orth(W)

  group.means <- group.means %*% W
  list(W=W, group.means=group.means, label=as.factor(clusters))
}

glda.pred <- function( train.x, train.y, test.x )
{
  y <- as.factor(train.y)
  n <- nrow(train.x)
  p <- ncol(train.x)
  c <- nlevels(y)
  if ( n < p )
  {
    model <- glda(train.x, y)
  } else
  {
    model <- lda(train.x, y)
  }
  pred <- knn(model$group.means, test.x%*%model$W, model$label)
  pred
}