NN.r 090624_2

# calculate FeedForward Network
# Online version

###########
# usage
# > source("NNn.r")
# > init(nin, nhidden, nout, fact, fout) -> some global variables are defined here
# > bprop(x, y, eta) -> learn by backpropagation
# > fprop(x) -> calculate feedforward propagation

# automatic usage
# > source("NNn.r")
# > x <- ... ( N*D dim. input matrix )
# > t <- ... ( N*K dim. output matrix )
# > out <- learn(x, t, hidden, times, eta, reset=FALSE)

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

###########
# functions
###########

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

###########
# arrange col-vectors of a matrix into a long column vector
###########

vec <- function(m){
  v <- NULL
  for(i in 1:dim(m)[2])
    v <- rbind(v, as.matrix(m[,i]))
  return ( as.matrix(v) )
}

###########
# rearrange vec into the original matrix
###########

mat <- function(v, ncol){
  nrow <- length(v)/ncol
  m <- NULL
  for(i in 1:ncol )
    m <- cbind(m, v[(((i-1)*nrow)+1):(i*nrow)])
  return (m)
}

###########
# "Declarations" of GLOBAL variables
# Once you init, the environment of R will be AFFECTED by the NN environment.
###########

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
f.out.d <- identity.d
Finv <- matrix(0)

###########
# initializer
###########

init <- function(nin, nhidden, nout, fact=sigmoid, fact.d=sigmoid.d, fout=identity, fout.d=identity.d){
  # 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
  f.out.d <- fout.d

  # 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

  Finv <<- diag( M*(D+1)+K*(M+1) )
}

###########
# set w1, w2
# renew D,M,K automatically.
# YOU SHOULD CHANGE "f." MANUALLY
###########

setw <- function(w1.in, w2.in){
  D <<- dim.in <<- (dim(w1.in)[2])-1
  M <<- dim.hidden <<- (dim(w1.in)[1])
  K <<- dim.out <<- (dim(w2.in)[1])

  w1 <<- w1.in
  w2 <<- w2.in

  Finv <<- diag( M*(D+1) + K*(M+1) )
}

###########
# calculate forward propergation
# x : D dim. input matrix
# note that x as well as z is always a vector, no more a matrix !
###########

fprop <- function( x ){
  # calculate hidden layers
  x.ex <- c(x, 1)
  a1 <- w1 %*% x.ex
  z <- f.act(a1)

  # calculate output layers
  z.ex <- c(z, 1)
  a2 <- w2 %*% z.ex
  y <- f.out(a2)

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

###########
# adaptive natural gradient descent - back propergation
# calculate natural gradients of w1 and d2 for a data n
# FOR sigmoid, identity ONLY! because each .d should be adjusted to their original functions only by rewriting the source below!
# data.in : input of training data
# data.out : output of training data
###########

aprop <- function(data.in, data.out, eta=0.01, rho=0.01){
  # FP
  calc <- fprop(data.in)

  # calculate output layers
  d2 <- calc$out - data.out
  e.d2 <- vec( d2 %o% c(calc$hidden, 1) )

  # calculate hidden layers
  # if the number of hidden nodes are only one, diag() would not work properly!
  if(M==1){
    d1 <- d2 %*% (w2[, -(M+1), drop=F]) %*% f.act.d(calc$hidden.a)
  } else {
    d1 <- d2 %*% (w2[, -(M+1), drop=F]) %*% diag( f.act.d(calc$hidden.a) )
  }
  e.d1 <- vec( as.vector(d1) %o% c(data.in, 1) )

  # renew parameters
  e.d <- rbind(e.d1, e.d2)   # the very order of e.d defined here determines the order of elements of Finv!!!
  Finv <<- ( Finv + rho * (Finv %*% e.d %*% t(e.d) %*% Finv ) / as.numeric( 1 + rho * t(e.d) %*% Finv %*% e.d ) )/( 1 - rho )
  v <- rbind( vec(w1), vec(w2) )
  v <- v - eta * Finv %*% e.d

  w1 <<- mat( v[1:(M*(D+1))], ncol=D+1)
  w2 <<- mat( v[(M*(D+1)+1):(K*(M+1))], ncol=M+1 )

  if(any(is.na(w2))){
    print("there's been founded some NAs in w2")
    return(list(d2, calc, error=FALSE))
  } else if(any(is.na(w1))) {
    print("there's been founded some NAs in w1")
    return(list(d1, data.in, error=FALSE))
  }

return ( list( w.hidden=w1, w.out=w2, error=sqrt(sum(d2^2)) ) )
}

###########
# back propergation
# calculate gradients of w1 and d2 for a data n
# FOR sigmoid, identity ONLY! because each .d should be adjusted to their original functions only by rewriting the source below!
# data.in : input of training data
# data.out : output of training data
###########

bprop <- function(data.in, data.out, eta=0.01){
  # FP
  calc <- fprop(data.in)

  # calculate output layers -> renew w2
  d2 <- calc$out - data.out
  w2 <<- w2 - eta * ( d2 %o% c(calc$hidden, 1) )

  # calculate hidden layers -> renew w1
  # if the number of hidden nodes are only one, diag() would not work properly!
  if(M==1){
    d1 <- d2 %*% (w2[, -(M+1), drop=F]) %*% f.act.d(calc$hidden.a)
  } else {
    d1 <- d2 %*% (w2[, -(M+1), drop=F]) %*% diag( f.act.d(calc$hidden.a) )
  }
    w1 <<- w1 - eta * ( as.vector(d1) %o% c(data.in, 1) )

  if(any(is.na(w2))){
    print("there's been founded some NAs in w2")
    return(list(d2, calc, error=FALSE))
  } else if(any(is.na(w1))) {
    print("there's been founded some NAs in w1")
    return(list(d1, data.in, error=FALSE))
  }

return ( list( w.hidden=w1, w.out=w2, error=sqrt(sum(d2^2)) ) )
}

###########
# learn from data
# both x and t must be data matrix (or, datum vector) in the form following:
#   x : N-row D-col Matrix
#   y : N-row K-col Matrix
# reset=FALSE : you use existing settings
# zeros=TRUE : initialize w1 and w2 with 0s
###########

learn <- function(x, t, hidden, times, eta, rho, reset=TRUE, zeros=FALSE){

  x <- as.matrix(x)
  t <- as.matrix(t)
  y <- NULL
  e <- NULL
  N <- dim(x)[1]

  if(reset) init(dim(x)[2], hidden, dim(t)[2])
  if(zeros){
    w1 <<- matrix(0, dim(w1)[1], dim(w1)[2])
    w2 <<- matrix(0, dim(w2)[1], dim(w2)[2])
  }

  for(i in 1:times){
    e.temp <- rep(0, N)
    for( n in 1:N ){
print(w1); print(w2); print(i);
      e.temp[n] <- aprop(x[n,],t[n,],eta,rho)$error
      if( !e.temp[n] ){
        print(i); print(n); print("error! learning stopped.")
        plot( c(e, sqrt(sum(e.temp^2))) )
        return( list( x=x, t=t, y=y, w1=w1, w2=w2, error=e ) )
      }
    }
y <- rbind(y, t(apply_fprop(x)))
    e <- c(e, sqrt(sum(e.temp^2))) # record errors on each cycle
  }
  plot(e)
#  return ( list( x=x, t=t, y=apply_fprop(x), w1=w1, w2=w2, error=e ) )
return ( list( x=x, t=t, y=y, w1=w1, w2=w2, error=e ) )
}

###########
# apply X to "fprop"
# x could be either a vector or a matrix.
# suppose X be the N*D dim matrix, as each data is arranged in row.
###########

apply_fprop <- function(x){
  if( !is.matrix(x) ) x <- t( as.matrix( x ) )

  temp <- NULL
  for( n in 1:(dim(x)[1]) ){
    result <- fprop( x[n,] )
    temp <- rbind(temp, result$out)
  }
  return(temp)
}

###########
# monitor hidden fire
# X is N*D dim matrix, expected which covers all valid inputs.
###########

monitor <- function(x){
  if( !is.matrix(x) ) x <- t( as.matrix( x ) )

  temp <- NULL
  for( n in 1:(dim(x)[1]) ){
    result <- fprop( x[n,] )
    temp <- rbind(temp, result$hidden.a)
  }
  return(temp)
}