# Neural Networks in R
# This code has been edited from 'http://euler.acadiau.ca/~hchipman/SpringSchool2006/nnet-intro.R'

library(MASS)
library(nnet)

#First example is XOR, the classes are 0,1
d<-data.frame(x=c(0,1,0,1),y=c(1,1,0,0))

#nnet defaults to using a least squares penalty function. 
t<-class.ind(c(1, 0, 0, 1))

mynet <-nnet(d[,1:2],t,size=2)
mynet.pred <- predict(mynet)

#outputs sums need not equal unity!
mynet.pred%*%c(1,1)

#What is the prediction for (0.5,0.5)?
predict(mynet,c(0.5,0.5))

#What function is the nnet actually evaluating?
#Generate a 100 by 100 lattice of points: 

m<-100
x<-seq(-0.1,1.1,length.out=m)
z<-data.frame(expand.grid(x=x,y=x))
pred<-predict(mynet,z)
#classes are 0, 1 so set contours at 0.5
plot(d[,1:2])
contour(x,x,matrix(max.col(pred),m,m),levels=c(1.5),
        add=T,d=F,lty=1,col=2)

#Second example: circle data
# simulate the data
p <- 2
N <- 200
set.seed(1)
x <- matrix(rnorm(N*p),ncol=p)
y <- as.numeric((x[,1]^2+x[,2]^2) > 1.4)
mydata <- data.frame(x=x,y=y)
plot(x,col=c('black','green')[y+1],pch=19,asp=1)

#It has one hidden layer, 3 hidden units, no weight decay, and is fit using "entropy",
#Note that we don't have to scale the inputs 
#(since they were generated from a distribution with mean 0 and standard deviation 1).

set.seed(3)

nn1 <- nnet(y~x.1+x.2,data=mydata,entropy=T,size=3,decay=0,maxit=2000,trace=T)

yhat <- as.numeric(predict(nn1,type='class'))
par(mfrow=c(1,2))
plot(x,pch=19,col=c('red','blue')[y+1],main='actual labels',asp=1)
plot(x,col=c('red','blue')[(yhat>0.5)+1],pch=19,main='predicted labels',asp=1)
table(actual=y,predicted=predict(nn1,type='class'))

# Re-run  the above with seed=4...
set.seed(4)
nn1 <- nnet(y~x.1+x.2,data=mydata,entropy=T,size=3,decay=0,maxit=2000,trace=T)
yhat <- as.numeric(predict(nn1,type='class'))
par(mfrow=c(1,2))
plot(x,pch=19,col=c('red','blue')[y+1],main='actual labels',asp=1)
plot(x,col=c('red','blue')[(yhat>0.5)+1],pch=19,main='predicted labels',asp=1)
table(actual=y,predicted=predict(nn1,type='class'))

# so the optimizer doesn't always find a good solution ...

# How many hidden units do we need?
par(mfrow=c(2,2))
for (i in 1:4){
  set.seed(3)
  nn1 <- nnet(y~x.1+x.2,data=mydata,entropy=T,size=i,decay=0,maxit=2000,trace=T)
  yhat <- as.numeric(predict(nn1,type='class'))
  plot(x,pch=20,col=c('black','green')[yhat+1])
  title(main=paste('nnet with',i,'hidden units'))
}

# let's look at what the function is actually doing.

set.seed(3)
nn1 <- nnet(y~x.1+x.2,data=mydata,entropy=T,size=3,decay=0,maxit=2000,trace=T)
x1grid <- seq(-3,3,l=200)
x2grid <- seq(-3,3,l=220)
xg <- expand.grid(x1grid,x2grid)
xg <- as.matrix(cbind(1,xg))
h1 <- xg%*%matrix(coef(nn1)[1:3],ncol=1)
h2 <- xg%*%matrix(coef(nn1)[4:6],ncol=1)
h3 <- xg%*%matrix(coef(nn1)[7:9],ncol=1)

par(mfrow=c(2,2))
contour(x1grid,x2grid,matrix(h1,200,220),levels=0)
contour(x1grid,x2grid,matrix(h2,200,220),levels=0,add=T)
contour(x1grid,x2grid,matrix(h3,200,220),levels=0,add=T)
title(main='boundaries from linear functions\n in hidden units')

sigmoid <- function(x){exp(x)/(1+exp(x))}
z <- coef(nn1)[10]+coef(nn1)[11]*sigmoid(h1)+coef(nn1)[12]*sigmoid(h2)+
  coef(nn1)[13]*sigmoid(h3)
contour(x1grid,x2grid,matrix(z,200,220))
title('sum of sigmoids \n of linear functions')
contour(x1grid,x2grid,matrix(sigmoid(z),200,220),levels=0.5)
title('sigmoid of sum of sigmoids \n of linear functions')
contour(x1grid,x2grid,matrix(sigmoid(z),200,220),levels=0.5)
points(x,pch=20,col=c('black','green')[y+1])
title('data values')

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

# EXERCISE:

# Work with a real-world data set: the 'glass' dataset:

# a dataset on fragments of glass collected in
# forensic work.  This data frame contains the following columns: RI
# refractive index; more precisely the refractive index is 1.518xxxx.
# The next 8 measurements are percentages by weight of oxides.  Na
# sodium.  Mg manganese.  Al aluminium.  Si silicon.  K potassium.  Ca
# calcium.  Ba barium.  Fe iron.  type The fragments were originally
# classed into seven types, one of which was absent in this dataset.

# The categories which occur are window float glass (WinF: 70), window non-float glass (WinNF: 76), vehicle window glass (Veh: 17), containers (Con: 13), tableware (Tabl: # 9) and vehicle headlamps (Head: 29).

# Hint code:

data (fgl)

mynet <- nnet (data=fgl,type ~ ., entropy=T,size=3,decay=0,maxit=2000,trace=T)

mynet.resubst <- predict (mynet, fgl, type='class')
mynet.resubst <- factor(mynet.resubst, levels=c("WinF", "WinNF", "Veh", "DCon", "Tabl", "Head"))

tab <- table(fgl$type, mynet.resubst)
mynet.resubst.error <- 1-sum(tab[row(tab)==col(tab)])/sum(tab)

# Aren'you satisfied with the result? Try the following:
# 1. scale() the data
# 2. increase the size (#hidden neurons)

# Design a validation protocol. Consider the possibility of performing 10-fold cross validation error estimation. Do we need a separate test set?

# needs the method 'cv10mlp{dprep}' (only for classification)
# suggestion: do it yourself! and take care of stratification!

# How do you select the number of hidden neurons ('size')
# Consider performing feature selection