NN.r 2

# calculate FeedForward Network

# usage
# > source("NN.r")
# > init(nin, nhidden, nout, fact, fout) -> some global variables are defined here
# > fprop(x) -> calculate feedforward propagation

# names below are reserved
# (---variables---)
# dim.in, dim.hidden, dim.out, D, M, K, w1, w2, f.act, f.act.d, f.out
# (---functions---)
# init, fprop, bprop, sigmoid, identity

# functions
sigmoid <- function(x) 1/(1+exp(-x))
sigmoid.d <- function(x) sigmoid(x)*( 1 - sigmoid(x) )
identity <- function(x) x

D <- numeric(0)
M <- numeric(0)
K <- numeric(0)
dim.in <- numeric(0)
dim.hidden <- numeric(0)
dim.out <- numeric(0)
w1 <- matrix(0)
w2 <- matrix(0)
f.act <- sigmoid
f.act.d <- sigmoid.d
f.out <- identity

# initializer
init <- function(nin, nhidden, nout, fact=sigmoid, fact.d=sigmoid.d, fout=identity){
# dimensions of input, hidden and output layers
D <<- dim.in <<- nin
M <<- dim.hidden <<- nhidden
K <<- dim.out <<- nout

f.act <<- fact
f.act.d <<- fact.d
f.out <<- fout

# initialize weights conforming to the PRML
w1 <<- matrix( runif( M*(D+1) ), M, D+1 ) # w1[,D+1] is bias
w2 <<- matrix( runif( K*(M+1) ), K, M+1 ) # w2[,M+1] is bias
}

# calculate forward propergation
# x : D*N dim. input matrix
fprop <- function( x ){
# calculate hidden layers
x <- as.matrix(x)   # x to be a 'column' vector
x.ex <- rbind( x, rep(1, dim(x)[2]) ) # bind (D+1)-th row to x
z <- f.act(w1 %*% x.ex)

# calculate output layers
z <- as.matrix(z)
z.ex <- rbind( z, rep(1, dim(z)[2]) )
y <- f.out(w2 %*% z.ex)

return( list( out=as.vector(y), hidden=as.vector(z) ) )
}

# back propergation
# calculate gradients of w1 and d2 for a data n
# data.in : input of training data
# data.out : output of training data
bprop <- function(data.in, data.out, eta=0.1){
# FP
old <- fprop(data.in)

# calculate output layers -> renew w2
d2 <- old$out - data.out
w2 <<- w2 - eta * ( as.vector(d2) %o% as.vector(old$hidden) )

# calculate hidden layers -> renew w1
d1 <- diag( f.act.d(old$out) ) %*% ( t(w2) %*% d2 )
w1 <<- w1 - eta * ( as.vector(d1) %o% as.vector(data.in) )

return(list(as.vector(d1), as.vector(d2)))
}

###########
# debug program
init(1,4,1)
x <- seq(0,3.14, by=0.1)
y <- sin(x)
for(i in 1:100){
for(n in 1:length(x))
  bprop(x[i],y[i])
}

p.y <- numeric(0)
for(p in x){
p.y <- c(p.y, fprop(p))
}

plot(x,p.y)