NN.r 完成版

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

###########
# "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

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

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

###########
# 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))
  } else if(any(is.na(w1))) {
    print("there's been founded some NAs in w1")
    return(list(d1, data.in))
  }

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
# if you set "reset=FALSE" you use existing settings
learn <- 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(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); 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)
}

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