source("pca.R")
source("rmvnorm.R")
source("mlUtilities.R")
source("fisherRRSVD.R")

n <- 200
X <- rmvnorm(matrix(rnorm(n*2),n,2), 0.3*matrix(c(2,1.8,1.8,2),2,2), c(2,3))
X <- rbind(X,rmvnorm(matrix(rnorm(n*2),n,2), 0.1*matrix(c(2,-1.3,-1.3,4),2,2), c(3,6)))

pr <- par(mfrow=c(1,2))

plot(X,col="red",pch=19,xlab=expression(x[1]),ylab=expression(x[2]),asp=1)

r <- pca(X)

means <- colMeans(X)
points(means[1],means[2],pch=18,cex=2,col="black")
Xmz <- X - matrix(rep(means,nrow(X)),nrow(X),ncol=2,byrow=TRUE)

drawVecs <- function (start,vs,labels=c(expression(v[1]),expression(v[2])),al=3) {
  arrows(start[1],start[2],start[1]+al*vs[1,],start[2]+al*vs[2,],lty=1,length=0.1,lwd=2)
  text(start[1]+(al+0.4)*vs[1,],start[2]+(al+0.4)*vs[2,],labels)
}
drawVecs(means,r$V,c(expression(e[1]),expression(e[2])),al=2)

newX <- Xmz %*% r$V

plot(newX,col="red",pch=19,xlab=expression(e[1]),ylab=expression(e[2]),asp=1)

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

### Read mpg data from UCI
mpg <-read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data")
mpgnames <- c("mpg","cylinders","displacement","horsepower",
              "weight","acceleration","year","origin","name")
colnames(mpg) <- mpgnames

### Remove all samples that have at least one "?"
keepRows <- apply(mpg != "?", 1, all)
mpg <- mpg[keepRows,]


### Assemble X and T matrices
X <- mpg[, 2:8]
X <- apply(X,2,as.numeric)
sf <- makeStandardizeF(X)
X <- sf(X)

pcaResult <- pca(X)
V <- pcaResult$V
rownames(V) <- colnames(mpg)[2:8]
colnames(V) <- paste("V",1:7,sep="")

x11()
pr <- par(mfrow=c(1,2))
plot(pcaResult$values,type="b",xlab="Indexs",ylab="Eigenvalue")

n<-nrow(V)
Vn <- matrix(t(V[7:1,]))
symbols(expand.grid(1:n,1:n),squares=abs(Vn), bg=ifelse(matrix(Vn)<0, "red", "green"), inches=FALSE,xaxt="n",yaxt="n",bty="n",xlab="",ylab="")
print(round(V,3))

######################################################################
### zip code data

drawSample <- function(row,X,T) {
  drawMatrix(X[row,])
  title(paste("Digit ",T[row]))
}

drawMatrix <- function(x) {
  image(matrix(x,16,16)[,16:1],xaxt="n",yaxt="n",xlab="",ylab="")
}


if (TRUE) {
  train <- as.matrix(read.table("zip.train"))
  X <- train[,-1]
  T <- train[,1,drop=FALSE]
}

pcaResult <- pca(X)

x11()
plot(pcaResult$values,type="l",ylab="Eigenvalues")

x11()
par(mfrow=c(3,3),mar=c(2,1,1,1))

for (i in 1:9) {
  drawMatrix(pcaResult$V[,i])
  title(paste("Eigenvector",i))
}

fisherResult <- fisherRRSVD(X,T,0)

x11()
par(mfrow=c(3,3),mar=c(2,1,1,1))

for (i in 1:9) {
  drawMatrix(fisherResult[,i])
  title(paste("Fisher vectors",i))
}

Xc <- X - matrix(colMeans(X),nrow(X),ncol(X),byrow=TRUE)

x11()
par(mfrow=c(1,2),mar=c(4,4,1,1))
plot(Xc %*% pcaResult$V[,1:2],pch=paste(T),col=T,type="p",xlab=expression(pca[1]),ylab=expression(pca[2]))
title("PCA")

plot(Xc %*% fisherResult[,1:2],pch=paste(T),col=T,type="p",xlab=expression(fish[1]),ylab=expression(fish[2]))
title("Fisher")