library("rgl")   ## for 3D graphics
library("e1071") ## for SVM

## load data from authors' website
# load(url(paste0("http://statweb.stanford.edu/~tibs/",
#                 "ElemStatLearn/datasets/ESL.mixture.rda")))
#dat  <- ESL.mixture

## load data from ESL R package
library("ElemStatLearn")
data(mixture.example)
dat <- mixture.example

datf <- data.frame(y=as.factor(dat$y), 
                   x1=dat$x[,1],
                   x2=dat$x[,2])

## fit SVM with radial kernel
## cost - the 'C' penalty
## gamma - parameter of the radial kernel function
## cost and gamma set to match HTF Fig 12.5
fit.svm <- svm(y ~ x1 + x2, kernel = "radial",
               cost = 1, gamma=1, data=datf)

## compute decision surface - f(x) 
dsf.svm <- predict(fit.svm, decision.values=TRUE,
                   newdata = data.frame(x1=dat$xnew[,1],
                                        x2=dat$xnew[,2]))
dsf.svm <- as.numeric(attr(dsf.svm, "decision.values"))
dat$dsf.svm <- dsf.svm
dat$dsf.svm <- with(dat, matrix(dsf.svm, length(px1), length(px2)))

## classification boundary - f(x) = 0
dat$cls.svm <- with(dat, contourLines(px1, px2, dsf.svm, levels=0))
## margins - f(x) = -1 and f(x) = 1
## because we set ||\beta|| = 1/M, margins occur at -1 and 1
dat$clsneg.svm <- with(dat, contourLines(px1, px2, dsf.svm, levels=-1))
dat$clspos.svm <- with(dat, contourLines(px1, px2, dsf.svm, levels=1))

## create 3D graphic, rotate to view 2D x1/x2 projection
par3d(FOV=1,userMatrix=diag(4))
plot3d(dat$xnew[,1], dat$xnew[,2], dat$prob, type="n",
       xlab="x1", ylab="x2", zlab="", zlim=range(dsf.svm),
       axes=FALSE, box=TRUE, aspect=1)

## plot points and bounding box; color margin support vectors black;
## fit.svm$coefs are estimates of 'alpha' multiplied by traning label (-1 or 1);
## the observations for which 0 < alpha < cost are support vectors
## that occur on the margin;
idx.svm <- fit.svm$index[which(abs(fit.svm$coefs) < 1)]
idx.svm <- 1:nrow(dat$x) %in% idx.svm
x1r <- range(dat$px1)
x2r <- range(dat$px2)
pts <- plot3d(dat$x[,1], dat$x[,2], max(dsf.svm),
              type="s", add=TRUE, 
              radius=ifelse(idx.svm, 0.05, 0.025), 
              col=ifelse(idx.svm, "black",
                         ifelse(dat$y, "orange", "blue")))
lns <- lines3d(x1r[c(1,2,2,1,1)], x2r[c(1,1,2,2,1)], max(dsf.svm))

## draw Bayes (True) classification boundary
dat$probm <- with(dat, matrix(prob, length(px1), length(px2)))
dat$cls <- with(dat, contourLines(px1, px2, probm, levels=0.5))
pls <- lapply(dat$cls, function(p) lines3d(p$x, p$y, z=max(dsf.svm),
                                           col="purple"))

## draw SVM classifier boundary and margins
pls <- lapply(dat$cls.svm, function(p) lines3d(p$x, p$y, z=max(dsf.svm), lwd=2))
plspos <- lapply(dat$clspos.svm, function(p) lines3d(p$x, p$y, z=max(dsf.svm)))
plsneg <- lapply(dat$clsneg.svm, function(p) lines3d(p$x, p$y, z=max(dsf.svm)))

## plot SVM surface and decision and margin planes
sfc <- surface3d(dat$px1, dat$px2, dsf.svm, alpha=1.0,
                 color="gray", specular="gray")
qds <- quads3d(x1r[c(1,2,2,1)], x2r[c(1,1,2,2)], 0, alpha=0.4,
               color="gray", lit=FALSE)
qdspos <- quads3d(x1r[c(1,2,2,1)], x2r[c(1,1,2,2)], 1, alpha=0.4,
                  color="gray", lit=FALSE)
qdsneg <- quads3d(x1r[c(1,2,2,1)], x2r[c(1,1,2,2)], -1, alpha=0.4,
                  color="gray", lit=FALSE)