# code based on files " Code regarding the first part of section 5.6",
# "Code regarding section 4.8", and "Code regarding section 5.7"
# (copyright 2003, 2004, 2012 A.Azzalini and B.Scarpa)

library(class)
library(e1071)

# help functions
pred.square <- function(model, x1, x2, grid = 50, Z.flag = FALSE)
{
  u <- pretty(x1, n = 10)
  x <- seq(u[1], u[10], length = grid)
  u <- pretty(x2, n = 10)
  y <- seq(u[1], u[10], length = grid)
  nx <- length(x)
  ny <- length(y)
  xoy <- cbind(rep(x, ny), as.vector(matrix(y, nx, ny, byrow = TRUE)))
  X <- matrix(xoy, nx * ny, 2, byrow = FALSE)
  if(Z.flag) pred <- predict(model, newdata = data.frame(z1 = X[ , 1], z2 = X[ , 2]))
  else  pred <- predict(model, newdata = data.frame(x1 = X[ , 1], x2 = X[ , 2]))
  list(pred = pred, x = x, y = y, X = X)
}


dataset <- read.table("http://azzalini.stat.unipd.it/Book-DM/classes.dat", header = TRUE, nrows = 200)
gr <- (2 - dataset[ , 3]) # groups are indicated by 0 and 1
K <- 2
n <- nrow(dataset)
attach(dataset)
n.grid <- 800
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]), ylab = expression(italic(x)[2]))
for(k in 1:K)
{
  g <- (group == k)
  points(dataset[g, 1:2], pch = 15 + k, col = 1 + k, cex = 1.5)
}

s3 <- svm(factor(group) ~ x1 + x2, data = dataset, kernel = "polynomial", degree = 3, cost = 1)
# figure 5.26a
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]), 
     ylab = expression(italic(x)[2]))
for(k in 1:K)
{
  g <- (group == k)
  points(dataset[g, 1:2], pch = 15 + k, col = 1 + k, cex = 1.5)
}
p <- pred.square(s3, x1, x2, grid = n.grid)
pred<- matrix(as.numeric(p$pred), n.grid, n.grid)
contour(p$x, p$y, pred, add = TRUE, levels = 1.5, lty = 2, drawlabels = FALSE)

#--------
s2 <- svm(factor(group) ~ x1 + x2, data = dataset, kernel = "radial", cost = 1)
# figure 5.26b
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]), ylab = expression(italic(x)[2]))
for(k in 1:K)
{
  g <- (group == k)
  points(dataset[g, 1:2], pch = 15 + k, col = 1 + k, cex = 1.5)
}
p <- pred.square(s2, x1, x2, grid = n.grid)
pred <- matrix(as.numeric(p$pred), n.grid, n.grid)
contour(p$x, p$y, pred, add = TRUE, levels = 1.5, lty = 2, drawlabels = FALSE)

s0P <- tune.svm(factor(group) ~ x1 + x2, data = dataset, kernel = "polynomial",
                cost = exp(seq(-1, 6, length = 25)))
# figure 5.26c
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]), ylab = expression(italic(x)[2]))
for(k in 1:K)
{
  g <- (group == k)
  points(dataset[g, 1:2], pch = 15 + k, col = 1 + k, cex = 1.5)
}
p <- pred.square(s0P$best.model, x1, x2, grid = n.grid)
pred <- matrix(as.numeric(p$pred), n.grid, n.grid)
contour(p$x, p$y, pred, add = TRUE, levels = 1.5, lty = 2, drawlabels = FALSE)

s0R <- tune.svm(factor(group) ~ x1 + x2, data = dataset, kernel = "radial",
                cost = exp(seq(-2, 4, length = 25)))
# figure 5.26d
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]), ylab = expression(italic(x)[2]))
for(k in 1:K)
{
  g <- (group == k)
  points(dataset[g, 1:2], pch = 15 + k, col = 1 + k, cex = 1.5)
}
p <- pred.square(s0R$best.model, x1, x2, grid = n.grid)
pred <- matrix(as.numeric(p$pred), n.grid, n.grid)
contour(p$x, p$y, pred, add = TRUE, levels = 1.5, lty = 2, drawlabels = FALSE)