######################################################################
### by Chuck Anderson for CS545, Fall 2009
### http://www.cs.colostate.edu/~anderson/cs545
######################################################################

source("nn.R")

### 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,]

### Randomly pick 80% of samples for training partition
randorder <- sample(nrow(mpg))
nTrain <- round(nrow(mpg)*0.8)
trainRows <- randorder[1:nTrain]

### Assemble X and T matrices
X <- mpg[trainRows, c(2,3,5,6,7)]
T <- mpg[trainRows, c(1,4), drop=FALSE]

### Convert all values to numeric
X <- apply(X,2,as.numeric)
T <- apply(T,2,as.numeric)

### Apply model w to test data
Xtest <- mpg[-trainRows, c(2,3,5,6,7)]
Ttest <- mpg[-trainRows, c(1,4), drop=FALSE]
Xtest <- apply(Xtest,2,as.numeric)
Ttest <- apply(Ttest,2,as.numeric)

rmse <- function(y,t) {
  sqrt(mean((y-t)^2))
}

result <- NULL
for (lam in c(0,0.001,0.01,0.1,1,10,20)) {
  for (nh in c(1,2,5,10,20,50)) {
    for (reps in 1:10) {
      nnet <- makeNN(X,T,nh,lambda=lam,nIter=10000)
      p <- useNN(nnet,Xtest)
      result <- rbind(result,
                      c(lam,nh,rmse(useNN(nnet,X),T),rmse(useNN(nnet,Xtest),Ttest)))
    }
  }
  means <- combine(result,1:2,3:4)
  print(means)
}

pr <- layout(matrix(c(1,1,2,3),2,2,byrow=TRUE))
lams <- sort(unique(means[,1]))
nhs <- sort(unique(means[,2]))
matresult <- matrix(means[,4],length(nhs),length(lams))
matplot(nhs, matresult,type="b",pch=paste(1:length(lams)),lty=1,xlab="Number of Nidden Units",ylab="Test RMSE",col=1:length(lams))
legend("topright",paste("a=",round(lams,3)),lty=1,col=1:length(lams),pch=paste(1:length(lams)))
persp(nhs,lams,matresult,xlab="Number of Hidden Units",ylab="Lambda",zlab="Test RMSE",
      theta=50,shade=0.8)
contour(nhs,lams,matresult,xlab="Number of Hidden Units",ylab="Lambda")

par(pr)

save(result,means,file="mpg.Rdata")