###########################################################################
## Copyright (c) 2002 Roger D. Peng
## 
## All rights reserved.
## 
## Permission is hereby granted, free of charge, to any person
## obtaining a copy of this software and associated documentation
## files (the "Software"), to deal in the Software without
## restriction, including without limitation the rights to use, copy,
## modify, merge, publish, distribute, and/or sell copies of the
## Software, and to permit persons to whom the Software is furnished
## to do so, provided that the above copyright notice(s) and this
## permission notice appear in all copies of the Software and that
## both the above copyright notice(s) and this permission notice
## appear in supporting documentation.
## 
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
## EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
## MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
## NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
## COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
## ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
## DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
## WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
## ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
## OF THIS SOFTWARE.
## 
## Except as contained in this notice, the name of a copyright holder
## shall not be used in advertising or otherwise to promote the sale,
## use or other dealings in this Software without prior written
## authorization of the copyright holder.
## 
###########################################################################

groupmatrix <- function(g) {
  h <- unique(g);
  ifelse(outer(g, h, "=="), 1, 0);
}


groupmean <- function(x, g) {
  h <- groupmatrix(g);
  d <- crossprod(h);
  dinv <- solve(d);
  dinv %*% crossprod(h, as.matrix(x));
}


## Compute the between groups covariance matrix

betweencov <- function(x, g) {
  h <- groupmatrix(g);
  dinv <- solve(crossprod(h));
  t(as.matrix(x)) %*% (h %*% (dinv %*% crossprod(h, as.matrix(x))));
}

betweencor <- function(x, g) {
  h <- groupmatrix(g);
  dinv <- solve(crossprod(h));
  t(as.matrix(x)) %*% (h %*% (dinv %*% crossprod(h, as.matrix(x))));
}




## Generalized eigenvalues and eigenvectors for matricies A and B.  Returns a 
## vector of eigenvalues and matrix of eigenvectors.

geneigen <- function(a, b) {
  s <- t(chol(b));
  tmat <- solve(s);
  m <- (tmat %*% a) %*% t(tmat);
  sv.decomp <- svd(m);
  eval <- sv.decomp$d^2;
  evecm <- sv.decomp$u;
  genevec <- crossprod(tmat, evecm);
  eigenvectors <- genevec;
  eigenvalues <- eval;

  list(values = eigenvalues, vectors = eigenvectors);
}


## Do a canonical discriminant analysis on data matrix DATA and group vector G.
## Returns a matrix of the canonical vectors, the generalized eigenvalues and a
## matrix of the generalized eigenvectors.

cda <- function(data, g, scale = F, do.print = F) {
  data <- scale(data, center=T, scale=scale);

  ct <- crossprod(data);
  cb <- betweencov(data, g);

  ngrp <- length(unique(g));
  geneig <- geneigen(cb, ct);
  canonicalvec <- data %*% geneig$vectors;

  if(do.print) {
    print(round(geneig$values[1:(ngrp - 1)], 4));
  }
  list(cvec = canonicalvec, values = geneig$values[1:(ngrp - 1)], 
       vectors = geneig$vectors, group = g);    
}