NN.r 090624

# 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, 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
identity.d <- function(x) return(1)

# arrange col-vectors 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 (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
F1 <- matrix(0)
F2 <- 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

  F1 <<- diag( M*(D+1) )
  F2 <<- diag( 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
}

###########
# 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
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)) ) )
}

###########
# 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, 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 ){
      e.temp[n] <- bprop(x[n,],t[n,],eta)$error
      if( !e.temp[n] ){
        print(i); print(n); print("error! learning stopped.")
        return
      }
    }
    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 ) )
}

###########
# !!! This method BADLY learn !!!
# learn from data, point by point.
# 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
# if you set "reset=FALSE" you use existing settings
learn_online <- function(x, t, hidden, times, eta, reset=TRUE){
  x <- as.matrix(x)
  t <- as.matrix(t)
  y <- NULL
  e <- NULL

  if(reset) init(dim(x)[2], hidden, dim(t)[2])

  for( n in 1:(dim(x)[1]) ){
    for(i in 1:times){
      e.out <- bprop(x[n,],t[n,],eta)
      if( !is.numeric(e.out$error) ){
        print(i); print(n); print("error! learning stopped.")
        return
      }
      e <- c(e, e.out$error)
    }
  }
  plot(e)
  return ( list( x=x, t=t, y=apply_fprop(x), 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)
}

###########
# debug program 1
# usage:
# > source("NN.r")
# > out <- test(times=2000, hidden=100, eta=0.01)
# > n <- 1:(dim(out$y)[1])    # means the number of loop counts
# > persp(n, out$x, out$y, theta=30, phi=30)

test <- function(times, hidden, eta){

  init(1,hidden,1)
  x <- seq(0, 3.14, by=0.1)
  t <- sin(x)
  e <- NULL
  y <- NULL

  for(i in 1:times){
    for(n in 1:length(x)){
      e.out <- bprop(x[n],t[n],eta)
      if(!is.numeric(e.out$error)){
        print(i,n)
        return
      }
      e <- c(e, e.out$error)
    }
  y <- rbind(y, apply_fprop(x))
  }

  plot(e)
  return (list(error=e, x=x, t=t, y=y))
}

###########
# debug program 2
# for simpler work (not return )

stest <- function(times, hidden, eta){

  init(1,hidden,1)
  x <- seq(0, 3.14, by=0.1)
  t <- sin(x)
  e <- NULL
  y <- NULL

  for(i in 1:times){
    for(n in 1:length(x)){
      e.out <- bprop(x[n],t[n],eta)
      if(!is.numeric(e.out$error)){
        print(i,n)
        return
      }
      e <- c(e, e.out$error)
    }
  }
  plot(e)
  return (list(error=e, x=x, t=t, y=apply_fprop(x)))
}

###########
# debug program 3
# BETE NOIRE XOR
#
# usage:
# > source("NNn.r")
# > out <- xor(times=2000, hidden=4, eta=0.1)
# > x <- y <- seq(-3,3,by=0.1)
# > z <- xor_view(x,y)

xor <- function(times, hidden, eta){
  init(2,hidden,1)
  x <- rbind( c(1,1), c(1,0), c(0,1), c(0,0) )
  t <- as.matrix(c(0, 1, 1, 0))
  e <- NULL
  y <- NULL

  for(i in 1:times){
    for( n in 1:(dim(x)[1]) ){
      e.out <- bprop(x[n,],t[n,],eta)
      if(!is.numeric(e.out$error)){
        print(i,n)
        return
      }
      e <- c(e, e.out$error)
    }
  }
  plot(e)
  return (list(error=e, x=x, t=t, y=apply_fprop(x)))
}

# for 3d view
# input x, y as vectors and show 3dgraph
xor_view <- function(x, y){
  z <- matrix( numeric(0), length(x), length(y))
  for(i in 1:length(x)){
    for(j in 1:length(y)){
      z[i,j] <- apply_fprop( c(x[i],y[j]))
    }
  }
  persp(x, y, z, theta=30, phi=30)
  return(z)
}