# Function gradPerc: Learn Single layer perceptron
# Input
# x - features (input)
# y - labels (output)
# g - learning rate
# it - maximum number of iterations
# Output
# w - perceptron weights
gradPerc <- function(x, y, g, it) {

  # transform data into appropriate format
  data <- rbind(1, x)

  # random initialisation of w = (theta, w1,..., wm) with uniform distribution
  ###### CODE HERE ######## -> Question 1
  w    <-  runif(nrow(data))

    # initialisation of vector for classification error
    err    <- numeric(it + 1)
  err[1] <- error(w, data, y)
  i      <- 1

  # run until all examples are classified correctly or until the maximum number of iterations
  while(err[i] > 0 & i <= it) {

    i <- i + 1

    # update the perceptron weights with gradient (batch) method and learning rate g
    w <- w + g * gradient(w, data, y)

    # normalise the weights
    if(w[1] != 0) {
      w <- w / abs(w[1])
    }

    # compute error in current iteration
    err[i] <- error(w, data, y)
  }

  # plot error over number of iterations
  plot(0:(i-1), err[1:i], main = "Errors", xlab = "iteration", ylab = "error",
       type = "l")

  return(w)
}

# Function err: compute number of classification errors
# Input
# w - perceptron weights
# x - features (input)
# y - labels (output)
# Output
# number of classification errors
error <- function(w, x, y) {
  # identify all incorrectly classified datapoints
  ###### CODE HERE ######## -> Question 2
  d = y * (t(w) %*% x)
  ind = d < 0

  # return number of classification errors
  ###### CODE HERE ######## -> Question 2
  return(sum(ind))
}

# Function gradient: compute gradient
# Input
# w - perceptron weights
# x - features (input)
# y - labels (output)
# Output
# gradient
gradient <- function(w, x, y) {
  # identify all incorrectly classified datapoints
  ###### CODE HERE ######## -> Question 3
  d = y * (t(w) %*% x)
  ind = d < 0

  # return gradient
  if(sum(ind) == 1) {
    ###### CODE HERE ######## -> Question 3
    return(y[ind] * x[,ind])
  } else {
    ###### CODE HERE ######## -> Question 3
    return ( rowSums(x[,ind] %*% diag(y[ind])))

  }
}

# read data from file
data1 <- read.table("data1.txt", header = TRUE, sep = "\t")
data2 <- read.table("data2.txt", header = TRUE, sep = "\t")

# Function learn_p_plot_data: Learn Single layer perceptron and plot data
# Input
# x - features (input)
# y - labels (output)
# g - learning rate
# it - maximum number of iterations
# Output
# w - perceptron weights
learn_p_plot_data <- function(x, y, g = 0.01, it = 100) {
  # learn perceptron
  w     <- gradPerc(t(x), y, 0.1, 100)

  # The following code adds a line to the scatterplot of the form y = a + b * x
  # Calculate a and b such that the line represents the separation threshold
  # defined by the learned weights of the perceptron

  ###### CODE HERE ######## -> Question 4
  a <-   w[1] / w[3] # line intercept
    b <-  -w[2] / w[3]  # line slope

    # plot data and seperator line
    plot(x, col = ifelse(y == 1, "blue", "red"), main = "Perceptron")
  abline(a, b)

  return(w)
}

# Code for Question 5:

set.seed(12)
q5_1 <- learn_p_plot_data(data1[, 1:2], data1[, 3])
q5_1
q5_2 <- learn_p_plot_data(data2[, 1:2], data2[, 3])
q5_2

# Code for Question 6:

set.seed(154)
q6_1 <- learn_p_plot_data(data1[, 1:2], data1[, 3], g = 0.01, it = 100)
q6_2 <- learn_p_plot_data(data1[, 1:2], data1[, 3], g = 0.10, it = 100)
q6_3 <- learn_p_plot_data(data1[, 1:2], data1[, 3], g = 1.00, it = 100)

q6_4 <- learn_p_plot_data(data2[, 1:2], data2[, 3], g = 0.01, it = 100)
q6_5 <- learn_p_plot_data(data2[, 1:2], data2[, 3], g = 0.10, it = 100)
q6_6 <- learn_p_plot_data(data2[, 1:2], data2[, 3], g = 1.00, it = 100)