require('kernlab')
dataA <- data.frame(
  x1 = c(0,4),
  x2 = c(0,4),
  y = c(1,-1)
)

dataA$y <- as.factor(dataA$y)

dataB <- data.frame(
  x1 = c(2,0,1),
  x2 = c(0,0,1),
  y = c(1,-1,-1)
)
dataB$y <- as.factor(dataB$y)

dataC <- data.frame(
  x1 = c(2, 2,-2,-2, 2, 2,-2,-2, 1, 1,-1,-1),
  x2 = c(2,-2,-2, 2, 2,-2,-2, 2, 1,-1,-1, 1),
  y = c( 1, 1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1)
)
dataC$y <- as.factor(dataC$y)

plot(dataA$x1, dataA$x2)

plot(dataB$x1, dataB$x2)

plot(dataC$x1, dataC$x2)

library(e1071)
formula <- y ~ .
svmA <- svm(formula, dataA, kernel="linear")
svmB <- svm(formula, dataB, kernel="linear")
svmC <- svm(formula, dataC, kernel="linear")

#Vectores de soporte
vsA <- dataA[svmA$index,1:2]
vsB <- dataB[svmB$index,1:2]
vsC <- dataC[svmC$index,1:2]
vsA
vsB
vsC

#Kernel
kA = matrix(c(0,0,0,32),2,2)
kB = matrix(c(4,0,1,0,0,0,1,0,2),3,3)

# Vector de pesos normal al hiperplano (W)
wA <- crossprod(as.matrix(vsA), svmA$coefs)
wB <- crossprod(as.matrix(vsB), svmB$coefs)
wC <- crossprod(as.matrix(vsC), svmC$coefs)
wA
wB
wC

# Calcular ancho del canal
widthA = 2/(sum(sqrt((wA)^2)))
widthB = 2/(sum(sqrt((wB)^2)))
widthC = 2/(sum(sqrt((wC)^2)))
widthA
widthB

# Calcular vector B
bA <- -svmA$rho
bB <- -svmB$rho
bC <- -svmC$rho
bA
bB
bC

# Calcular la ecuación del hiperplano y de los planos de soporte positivo y negativo
paste(c("[",wA,"]' * x + [",bA,"] = 0"), collapse=" ")
paste(c("[",wA,"]' * x + [",bA,"] = 1"), collapse=" ")
paste(c("[",wA,"]' * x + [",bA,"] = -1"), collapse=" ")

paste(c("[",wB,"]' * x + [",bB,"] = 0"), collapse=" ")
paste(c("[",wB,"]' * x + [",bB,"] = 1"), collapse=" ")
paste(c("[",wB,"]' * x + [",bB,"] = -1"), collapse=" ")

paste(c("[",wC,"]' * x + [",bC,"] = 0"), collapse=" ")
paste(c("[",wC,"]' * x + [",bC,"] = 1"), collapse=" ")
paste(c("[",wC,"]' * x + [",bC,"] = -1"), collapse=" ")

# Determinar la clase a la que pertenece un punto dado
x = c(5,6)
if ((t(wA) %*% x + bA) >= 0){
  print("1")
}else if((t(wA) %*% x + bA) < 0){
  print("-1")
}

x = c(1,4)
if ((t(wA) %*% x + bA) >= 0){
  print("1")
}else if((t(wA) %*% x + bA) < 0){
  print("-1")
}

x = c(5,6)
if ((t(wB) %*% x + bB) >= 0){
  print("1")
}else if((t(wB) %*% x + bB) < 0){
  print("-1")
}

x = c(1,4)
if ((t(wB) %*% x + bB) >= 0){
  print("1")
}else if((t(wB) %*% x + bB) < 0){
  print("-1")
}

x = c(1,0)
if ((t(wC) %*% x + bC) >= 0){
  print("1")
}else if((t(wC) %*% x + bC) < 0){
  print("-1")
}

x = c(-1,0)
if ((t(wC) %*% x + bC) >= 0){
  print("1")
}else if((t(wC) %*% x + bC) < 0){
  print("-1")
}

#Dibujar los puntos y el hiperplano
plot(svmA, dataA)

plot(svmB, dataB)

plot(svmC, dataC)

kfunction <- function(x1 =0, x2=0){
  if(sqrt(x1^2+x2^2)>2){
    k <- function(x1,x2){
      c(4-x2+ abs(x1-x2),4-x1+ abs(x1-x2))
    }
  }
  else{
    k <- function(x1,x2){
      c(x1,x2)
    }
  }
  class(k) <- "kernel"
  k
}

x1 = matrix(cbind(dataC$x1,dataC$x2),2)
svp <- ksvm(x1,dataC$y,type="C-svc",C = 100, kernel=kfunction(1,0),scaled=c())

plot(c(min(x1[,1]), max(x1[,1])),c(min(x1[,2]), max(x1[,2])),type='n',xlab='x1',ylab='x2')
title(main='Linear Separable Features')
ymat <- ymatrix(svp)
points(x1[-SVindex(svp),1], x1[-SVindex(svp),2], pch = ifelse(ymat[-SVindex(svp)] < 0, 2, 1))
points(x1[SVindex(svp),1], x1[SVindex(svp),2], pch = ifelse(ymat[SVindex(svp)] < 0, 17, 16))

# Extract w and b from the model
w <- colSums(coef(svp)[[1]] * x1[SVindex(svp),])
b <- b(svp)

# Draw the lines
abline(b/w[2],-w[1]/w[2])
abline((b+1)/w[2],-w[1]/w[2],lty=2)
abline((b-1)/w[2],-w[1]/w[2],lty=2)

# Determinar la clase a la que pertenece un punto dado
x = c(4,5)
if ((t(w) %*% x + b) >= 0){
  print("1")
}else if((t(w) %*% x + b) < 0){
  print("-1")
}

svmIris <- svm(formula, iris, kernel="linear")

#Vectores de soporte
vsIris <- dataA[svmIris$index,1:2]
vsIris

# Vector de pesos normal al hiperplano (W)
wIris <- crossprod(as.matrix(vsIris), svmIris$coefs)
wIris

# Calcular ancho del canal
widthIris = 2/(sum(sqrt((wIris)^2)))
widthIris

# Calcular vector B
bIris <- -svmIris$rho
bIris

# Calcular la ecuación del hiperplano y de los planos de soporte positivo y negativo

paste(c("[",wIris,"]' * x + [",bIris,"] = 0"), collapse=" ")
paste(c("[",wIris,"]' * x + [",bIris,"] = 1"), collapse=" ")
paste(c("[",wIris,"]' * x + [",bIris,"] = -1"), collapse=" ")