###
###
###  Newral Network Demo 1 - Spiral Response
###
###
####################################################
### https://selbydavid.com/2018/01/09/neural-network/


sigmoid <- function(x) 1 / (1 + exp(-x))

feedforward <- function(x, w1, w2) {
  z1 <- cbind(1, x) %*% w1
  h <- sigmoid(z1)
  z2 <- cbind(1, h) %*% w2
  list(output = sigmoid(z2), h = h)
}

backpropagate <- function(x, y, y_hat, w1, w2, h, learn_rate) {
  dw2 <- t(cbind(1, h)) %*% (y_hat - y)
  dh  <- (y_hat - y) %*% t(w2[-1, , drop = FALSE])
  dw1 <- t(cbind(1, x)) %*% (h * (1 - h) * dh)

  w1 <- w1 - learn_rate * dw1
  w2 <- w2 - learn_rate * dw2

  list(w1 = w1, w2 = w2)
}

train <- function(x, y, hidden = 5, learn_rate = 1e-2, iterations = 1e4) {
  d <- ncol(x) + 1
  w1 <- matrix(rnorm(d * hidden), d, hidden)
  w2 <- as.matrix(rnorm(hidden + 1))
  for (i in 1:iterations) {
    ff <- feedforward(x, w1, w2)
    bp <- backpropagate(x, y,
                        y_hat = ff$output,
                        w1, w2,
                        h = ff$h,
                        learn_rate = learn_rate)
    w1 <- bp$w1; w2 <- bp$w2
  }
  list(output = ff$output, w1 = w1, w2 = w2)
}



##--- Simulate Spiral Data Set
two_spirals <- function(N = 200,
                        radians = 3*pi,
                        theta0 = pi/2,
                        labels = 0:1) {
  N1 <- floor(N / 2)
  N2 <- N - N1

  theta <- theta0 + runif(N1) * radians
  spiral1 <- cbind(-theta * cos(theta) + runif(N1),
                   theta * sin(theta) + runif(N1))
  spiral2 <- cbind(theta * cos(theta) + runif(N2),
                   -theta * sin(theta) + runif(N2))

  points <- rbind(spiral1, spiral2)
  classes <- c(rep(0, N1), rep(1, N2))

  data.frame(x1 = points[, 1],
             x2 = points[, 2],
             class = factor(classes, labels = labels))
}
set.seed(42)
hotdogs <- two_spirals(labels = c('not hot dog', 'hot dog'))

library(ggplot2)
theme_set(theme_classic())
ggplot(hotdogs) +
  aes(x1, x2, colour = class) +
  geom_point() + labs(x = expression(x[1]), y = expression(x[2]))




##--- Training NN
x <- data.matrix(hotdogs[, c('x1', 'x2')])
y <- hotdogs$class == 'hot dog'


 nnet5 <- train(x, y, hidden = 10, iterations = 1e5)
# nnet5 <- train(x, y, hidden = 30, iterations = 1e5)
# nnet5 <- train(x, y, hidden = 1, iterations = 1e5)

mean((nnet5$output > .5) == y)

grid <- expand.grid(x1 = seq(min(hotdogs$x1) - 1,
                             max(hotdogs$x1) + 1,
                             by = .25),
                    x2 = seq(min(hotdogs$x2) - 1,
                             max(hotdogs$x2) + 1,
                             by = .25))

ff_grid <- feedforward(x = data.matrix(grid[, c('x1', 'x2')]),
                       w1 = nnet5$w1,
                       w2 = nnet5$w2)
grid$class <- factor((ff_grid$output > .5) * 1,
                     labels = levels(hotdogs$class))

ggplot(hotdogs) + aes(x1, x2, colour = class) +
  geom_point(data = grid, size = .5) +
  geom_point() +
    labs(x = expression(x[1]), y = expression(x[2]))