require(rgl)   #Functions for interactive 3D-plotting
require(e1071) #Contains an implementation of the algorithm for estimation of support vector machines
require(MASS)  #Contains the dataset iris used
###require(kernlab) #Another implementation of support vector machines through the function ksvm

###Data load and setup 

data(iris)

spec <- data.frame(setosa=1:50,versicolor=51:100,virginica=101:150)

###Helper functions for 2D-plotting of the data

plot2 <- function(i,j,species=c("setosa","versicolor"))
{
  data1 <- iris[spec[species[1]][[1]],c(i,j)]
  data2 <- iris[spec[species[2]][[1]],c(i,j)]
    
  plot(data1,
       xlim=range(c(data1[,1],data2[,1])) + c(-1,1),
       ylim=range(c(data1[,2],data2[,2])) + c(-1,1),
       col="red",pch=19)
  
  points(data2,col="blue",pch=19)
  return(c(i,j))
}

addLine <- function(i,j,species=c("setosa","versicolor"))
{
  ind <- c(i,j)
  specInd <- (1:150)[iris[,5] %in% species]
  xData <- iris[specInd,ind]

  svmOb <- svm(x=xData,y=iris[specInd,5],
               type="C-classification",
               kernel="linear",
               cost=1000,scale=FALSE)
  
  betaHat <- t(svmOb$coefs) %*% as.matrix(iris[specInd[svmOb$ind],ind])
  beta0 <- svmOb$rho
  f0 <- function(x) {-betaHat[1]*x/betaHat[2] + beta0/betaHat[2]}

  lines(range(xData[,1]),sapply(range(xData[,1]),f0))
  points(xData[svmOb$ind,],pch=19,cex=2)
}

addKnn <- function(i,j,species=c("setosa","versicolor"),k=5){
  ind <- c(i,j)
  specInd <- (1:150)[iris[,5] %in% species]
  xData <- iris[specInd,ind]
  
  knnClass <- function(x,y){
    knn(xData,c(x,y),cl=iris[specInd,5],k=k)
  }

  xx <- seq(min(xData[,1])-1,max(xData[,1])+1,length.out=100)
  yy <- seq(min(xData[,2])-1,max(xData[,2])+1,length.out=100)
  knnResult <- matrix(nrow=length(xx),ncol=length(yy))

  for(i in 1:length(xx)){
    for(j in 1:length(yy)){
      knnResult[i,j] <- knnClass(xx[i],yy[j])
    }
  }
  
  contour(xx,yy,knnResult,levels=c(1,2,3),xlim=c(4,8),ylim=c(0,7),add=TRUE)
}
  


###Helper functions for 3D-plotting of the data

plot3 <- function(i,j,k,species=c("setosa","versicolor"))
{
  data1 <- iris[spec[species[1]][[1]],c(i,j,k)]
  data2 <- iris[spec[species[2]][[1]],c(i,j,k)]
  
  plot3d(data1,
         xlim=range(c(data1[,1],data2[,1])) +c(-1,1),
         ylim=range(c(data1[,2],data2[,2])) +c(-1,1),
         zlim=range(c(data1[,3],data2[,3])) +c(-1,1),
         type="s",col="red",radius=0.05)
  
  plot3d(data2,type="s",col="blue",radius=0.05,add=TRUE)

  return(c(i,j,k))
}

addPlane <-function(i,j,k,species=c("setosa","versicolor"))
{
  ind <- c(i,j,k)
  specInd <- (1:150)[iris[,5] %in% species]
  xData <- iris[specInd,ind]
  
  svmOb <- svm(x=xData,y=iris[specInd,5],
               type="C-classification",
               kernel="linear",
               scale=FALSE,cost=1000)
  
  betaHat <- t(svmOb$coefs) %*% as.matrix(xData[svmOb$ind,])
  beta0 <- svmOb$rho
  f0 <- function(x) {-betaHat[1]*x[1]/betaHat[3] - betaHat[2]*x[2]/betaHat[3] + beta0/betaHat[3]}

  x <- range(xData[,1]) 
  y <- range(xData[,2])
  z <- c(f0(c(x[1],y[1])),f0(c(x[1],y[2])),f0(c(x[2],y[2])),f0(c(x[2],y[1])))
  z0 <- range(xData[,3]) + c(-1,1)
  rangeOK <- (min(z) >= z0[1]) & (max(z) <= z0[2])
  
  planeCoords <- data.frame(x=c(x[1],x[1],x[2],x[2]),y=c(y,y[2],y[1]),z=z)
  quads3d(planeCoords,alpha=0.3)
  
  plot3d(xData[svmOb$ind[svmOb$ind <=50],],
         type="s",col="red",radius=0.1,add=TRUE)
  
  plot3d(xData[svmOb$ind[svmOb$ind > 50],],
         type="s",col="blue",radius=0.1,add=TRUE)

}

###
###Examples
###


plot2(1,2)
addLine(1,2)
plot2(1,2,c("versicolor","virginica"))
addKnn(1,2,c("versicolor","virginica"))
addLine(1,2,c("versicolor","virginica"))

plot3(1,2,3)
addPlane(1,2,3)


plot3(1,2,3,c("versicolor","virginica"))
addPlane(1,2,3,c("versicolor","virginica"))