library(MASS)
library(nnet)
library(class)
zip.train <- read.table("https://www.uio.no/studier/emner/matnat/math/STK2100/data/zip.train",
                        col.names=c("cl",paste("x",1:256,sep="")))
zip.test <- read.table("https://www.uio.no/studier/emner/matnat/math/STK2100/data/zip.test",
                       col.names=c("cl",paste("x",1:256,sep="")))
zip = rbind(zip.train,zip.test)
N.tr = nrow(zip.train)
N.te = nrow(zip.test)
#Normalizing covariates
zip[,-1] = scale(zip[,-1])
                       
par(mfrow=c(4,4),mar=c(2,2,2,2))
for(i in 1:16)
  {
  	foo = matrix(as.matrix(zip.train[i,-1]),ncol=16)  
  	image(foo[,16:1],main=zip.train$cl[i],col=gray((0:32)/32),xaxt="n",yaxt="n")
  }                     
zip$cl = as.factor(zip$cl)
zip.train = zip[1:N.tr,]
zip.test = zip[N.tr+1:N.te,]
zip.train$cl = as.factor(zip.train$cl)
zip.test$cl = as.factor(zip.test$cl)

#Neural network with one hidden layer
#Note: The following call is VERY time-consuming
decay = seq(0,4.5,length=10)
err.nnet = rep(NA,length(decay))
##PRERUN GAVE THE FOLLOWING ERROR RATES
err.nnet = 
 c(0.15097160, 0.09417040, 0.08719482, 0.08669656, 0.08221226, 
   0.08071749, 0.08719482, 0.08619831, 0.09118087, 0.07922272)
if(0)
{
for(i in 1:length(decay))
{
 m = 10
 zip.nnet = nnet(cl~.,data=zip[1:N.tr,],size=m,decay=decay[i],
   MaxNWts=10000,maxit=500)
 pred = predict(zip.nnet,zip[N.tr+1:N.te,],type="class")
 table(zip.test$cl,pred)
 err.nnet[i]=mean(zip.test$cl!=pred)
}
}
par(mfrow=c(1,1))
plot(decay,err.nnet)

#Deep learning with three hidden layers
library(RSNNS)
zipTargets <- decodeClassLabels(zip[,1])
zipTargets.train = zipTargets[1:N.tr,]
zipTargets.test = zipTargets[N.tr+1:N.te,]
maxit = c(100,200,300,400,500,600,700,800,900,1000)
err.dnet = rep(NA,length(maxit))
#PRERUN
err.dnet = c(0.1305431, 0.1335326, 0.1145989, 0.1175884, 0.1275536, 
             0.1190832, 0.1081216, 0.1285501, 0.1170902, 0.1320379)
if(0)
{
for(i in 1:length(maxit))
{
zip.dnet = mlp(as.matrix(zip.train[,-1]),zipTargets.train, size = c(6,10,6),
               learnFuncParams=c(0.3),maxit=maxit[i])
pred.dnet = apply(predict(zip.dnet,zip.test[,-1]),1,which.is.max)-1
#table(zip.test$cl,pred.dnet)
err.dnet[i] = mean(zip.test$cl!=pred.dnet)
}
}
plot(maxit,err.dnet)
zip.dnet = mlp(as.matrix(zip.train[,-1]),zipTargets.train, size = c(6,10,6),
               learnFuncParams=c(0.3),maxit=1000)
pred = predict(zip.dnet,zip.test[,-1])
pred.dnet = apply(pred,1,which.is.max)-1
#table(zip.test$cl,pred.dnet)
show(mean(zip.test$cl!=pred.dnet))

#Deep learning using first 10 principal component
#Mainly to show how to display the final network
zip.pca = prcomp(zip[,-1],retx=TRUE)$x
zip.dnet2 = mlp(as.matrix(zip.pca[1:N.tr,1:10]),zipTargets.train, size = c(6,10,6),
                learnFuncParams=c(0.3),
                maxit=500, inputsTest=zip.pca[N.tr+1:N.te,1:10], targetsTest=zipTargets.test)
pred.dnet2 = apply(zip.dnet2$fittedTestValues,1,which.max)-1
table(zip.test$cl,pred.dnet2)
mean(zip.test$cl!=pred.dnet2)
summary(zip.dnet2)

#import plot.nnet function from Github
library(devtools)
source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r')

#Make plot
par(mar=numeric(4),mfrow=c(1,1),family='serif')
plot.nnet(zip.dnet2)