# Paleta: clase 1 = morado, clase 2 = azul, clase 3 = naranja
pal <- c("#C45BC4", "#2BA8E0", "#E8A33D")

# Puntaje de un SVM binario: > 0  <=>  segunda clase (levels[2]).
# (e1071 devuelve el signo contrario: positivo = primera clase)
score_svm <- function(fit, newdata) {
  dv <- attr(predict(fit, newdata, decision.values = TRUE), "decision.values")
  -as.numeric(dv[, 1])
}

# Coeficientes de un SVM LINEAL (solo si scale = FALSE):  f(x) = beta0 + beta'x
# e1071 no los entrega directamente; se reconstruyen desde los vectores de soporte.
coef_svm <- function(fit) {
  stopifnot(fit$kernel == 0, !any(fit$scaled))
  beta <- -colSums(fit$coefs[, 1] * fit$SV)   # signo: positivo = levels[2]
  list(beta0 = fit$rho, beta = beta, margen = 1 / sqrt(sum(beta^2)))
}

# Frontera de decisión (línea continua), margen (punteada) y vectores de soporte (círculos)
plot_svm <- function(fit, data, titulo = NULL, margenes = TRUE, n = 150) {
  lev <- levels(data$y)
  rejilla <- expand_grid(
    x1 = seq(min(data$x1) - .4, max(data$x1) + .4, length.out = n),
    x2 = seq(min(data$x2) - .4, max(data$x2) + .4, length.out = n))
  rejilla$pred <- predict(fit, rejilla)
  binario <- length(lev) == 2
  if (binario) rejilla$score <- score_svm(fit, rejilla)
  data$sv <- seq_len(nrow(data)) %in% fit$index

  p <- ggplot() +
    geom_raster(data = rejilla, aes(x1, x2, fill = pred), alpha = .15) +
    geom_point(data = data, aes(x1, x2, colour = y), size = 2) +
    geom_point(data = filter(data, sv), aes(x1, x2), shape = 1, size = 4.5) +
    scale_fill_manual(values = pal, guide = "none") +
    scale_colour_manual(values = pal) +
    coord_cartesian(expand = FALSE) +
    labs(title = titulo, colour = "clase", x = expression(X[1]), y = expression(X[2]))
  if (binario) {
    p <- p + geom_contour(data = rejilla, aes(x1, x2, z = score),
                          breaks = 0, colour = "black", linewidth = .8)
    if (margenes)
      p <- p + geom_contour(data = rejilla, aes(x1, x2, z = score),
                            breaks = c(-1, 1), colour = "black",
                            linetype = "dashed", linewidth = .4)
  }
  p
}

# Matriz de confusión como mapa de calor
plot_conf <- function(pred, real, titulo = NULL) {
  tibble(predicho = pred, real = real) |>
    count(predicho, real) |>
    ggplot(aes(real, predicho, fill = n)) +
    geom_tile(colour = "white") + geom_text(aes(label = n), size = 6) +
    scale_fill_gradient(low = "white", high = pal[2], guide = "none") +
    labs(title = titulo, x = "clase real", y = "clase predicha")
}

# Curvas ROC con pROC (lista nombrada de objetos roc)
plot_roc <- function(lst, titulo = NULL) {
  nombres <- imap_chr(lst, \(r, nm) sprintf("%s (AUC = %.3f)", nm, as.numeric(auc(r))))
  ggroc(setNames(lst, nombres), legacy.axes = TRUE, linewidth = 1) +
    geom_abline(linetype = 3) + coord_equal() +
    labs(title = titulo, x = "Tasa de falsos positivos",
         y = "Tasa de verdaderos positivos", colour = NULL) +
    theme(legend.position = "bottom", legend.direction = "vertical")
}

# Distancia con signo y pie de la perpendicular
dist_h <- function(x0, beta0, beta) (beta0 + sum(beta * x0)) / sqrt(sum(beta^2))
pie_h  <- function(x0, beta0, beta) x0 - (beta0 + sum(beta * x0)) / sum(beta^2) * beta

b0 <- 1; b <- c(2, 3)          # el hiperplano de la Figura 9.1
P  <- c(1, 1)                  # un punto
Q  <- pie_h(P, b0, b)

c(distancia_con_signo = dist_h(P, b0, b),
  a_mano              = 6 / sqrt(13),                # |1 + 2 + 3| / sqrt(4 + 9)
  norma_de_P_menos_Q  = sqrt(sum((P - Q)^2)),
  Q_en_el_hiperplano  = b0 + sum(b * Q))             # debe ser 0

# Verificación "de cálculo": minimizar la distancia sobre la recta, a ciegas
sobre_H <- \(t) c(t, -(b0 + b[1] * t) / b[2])        # parametrización de H
optimize(\(t) sum((sobre_H(t) - P)^2), c(-5, 5))$objective |> sqrt()

# Simula dos nubes y repite con otra semilla hasta que sean separables
sim_separable <- function(seed = 1, n = 20, shift = c(2.5, -2.5)) {
  repeat {
    set.seed(seed)
    x <- matrix(rnorm(n * 2), ncol = 2)
    y <- rep(c(-1, 1), each = n / 2)
    x[y == 1, ] <- sweep(x[y == 1, ], 2, shift, "+")
    d <- tibble(x1 = x[, 1], x2 = x[, 2], y = factor(y))
    ajuste <- svm(y ~ ., data = d, kernel = "linear", cost = 1e5, scale = FALSE)
    if (all(fitted(ajuste) == d$y)) return(d)
    seed <- seed + 1
  }
}
d_sep <- sim_separable()
mm    <- svm(y ~ ., data = d_sep, kernel = "linear", cost = 1e5, scale = FALSE)
cf    <- coef_svm(mm)
cf

yy <- if_else(d_sep$y == levels(d_sep$y)[2], 1, -1)
X  <- as.matrix(d_sep[c("x1", "x2")])
c0 <- pie_h(colMeans(X), cf$beta0, cf$beta)          # centro de la "brecha", sobre H
rot <- \(th, v) c(cos(th) * v[1] - sin(th) * v[2], sin(th) * v[1] + cos(th) * v[2])

alt <- tibble(grados = seq(-60, 60, by = 5)) |>
  mutate(beta  = map(grados * pi / 180, rot, v = cf$beta),
         beta0 = map_dbl(beta, \(b) -sum(b * c0)),
         b1 = map_dbl(beta, 1), b2 = map_dbl(beta, 2),
         separa = map2_lgl(beta, beta0, \(b, b0) all(yy * (b0 + X %*% b) > 0)))

sep    <- filter(alt, separa)
lineas <- bind_rows(slice(sep, c(1, n())), filter(alt, grados == 0))

p_izq <- ggplot(d_sep, aes(x1, x2, colour = y)) +
  geom_abline(data = lineas, aes(intercept = -beta0 / b2, slope = -b1 / b2)) +
  geom_point(size = 2.5) + scale_colour_manual(values = pal) +
  coord_cartesian(xlim = range(d_sep$x1), ylim = range(d_sep$x2)) +
  labs(title = "Tres hiperplanos separadores")
p_der <- plot_svm(mm, d_sep, "Uno de ellos define la regla", margenes = FALSE)
p_izq + p_der + plot_layout(guides = "collect")

plot_svm(mm, d_sep, sprintf("Margen máximo M = %.3f; %d vectores de soporte", cf$margen, nrow(mm$SV)))

# Agregamos un punto de la clase 2 en el centro de la brecha
p_nuevo <- c0 + 0.1 * cf$beta / sqrt(sum(cf$beta^2))
d_nuevo <- bind_rows(d_sep, tibble(x1 = p_nuevo[1], x2 = p_nuevo[2],
                                   y = factor(levels(d_sep$y)[2], levels = levels(d_sep$y))))
mm2 <- svm(y ~ ., data = d_nuevo, kernel = "linear", cost = 1e5, scale = FALSE)

plot_svm(mm,  d_sep,   sprintf("Sin el punto: M = %.2f", cf$margen), margenes = FALSE) +
plot_svm(mm2, d_nuevo, sprintf("Con el punto: M = %.2f", coef_svm(mm2)$margen), margenes = FALSE) +
  plot_layout(guides = "collect")

set.seed(2)
d_sol <- tibble(x1 = rnorm(40), x2 = rnorm(40), y = factor(rep(c(-1, 1), each = 20))) |>
  mutate(x1 = x1 + 1.2 * (y == "1"), x2 = x2 - 1.2 * (y == "1"))

ggplot(d_sol, aes(x1, x2, colour = y)) + geom_point(size = 3) +
  scale_colour_manual(values = pal) + coord_fixed()

aj <- svm(y ~ ., data = d_sol, kernel = "linear", cost = 1, scale = FALSE)

tabla <- d_sol |>
  mutate(f        = score_svm(aj, d_sol),
         yy       = if_else(y == "1", 1, -1),
         holgura  = pmax(0, 1 - yy * f),
         ubicacion = case_when(
           holgura > 1    ~ "3. lado equivocado del hiperplano (ε > 1)",
           holgura > 1e-6 ~ "2. viola el margen (0 < ε ≤ 1)",
           TRUE           ~ "1. lado correcto del margen (ε = 0)"))

count(tabla, ubicacion)

c(suma_holguras = sum(tabla$holgura), n_vectores_soporte = length(aj$index))

costos <- c(0.01, 0.1, 1, 100)
graficas <- map(costos, \(k) {
  aj <- svm(y ~ ., data = d_sol, kernel = "linear", cost = k, scale = FALSE)
  plot_svm(aj, d_sol, sprintf("cost = %g  ·  %d vectores de soporte", k, nrow(aj$SV)))
})
wrap_plots(graficas, ncol = 2, guides = "collect")

set.seed(1)
d_nl <- bind_rows(
  tibble(x1 = rnorm(50, 0.5, 0.9), x2 = rnorm(50, -0.2, 0.9), y = "A"),
  tibble(x1 = rnorm(25, 2.8, 1),   x2 = rnorm(25, 2.8, 1),    y = "B"),
  tibble(x1 = rnorm(25, -3, 1),    x2 = rnorm(25, -3, 1),     y = "B")) |>
  mutate(y = factor(y))

svc_lin <- svm(y ~ ., data = d_nl, kernel = "linear", cost = 1, scale = FALSE)

p1 <- ggplot(d_nl, aes(x1, x2, colour = y)) + geom_point(size = 2) +
  scale_colour_manual(values = pal) + labs(title = "Datos")
p2 <- plot_svm(svc_lin, d_nl, "Clasificador de vectores de soporte (lineal)")
p1 + p2 + plot_layout(guides = "collect")

svc_cuad <- svm(y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2), data = d_nl,
                kernel = "linear", cost = 1)
plot_svm(svc_cuad, d_nl, "SVC lineal en (X1, X2, X1², X2², X1·X2)")

mean(predict(svc_cuad, d_nl) == d_nl$y)    # exactitud de entrenamiento

phi <- \(x) c(1, sqrt(2) * x[1], sqrt(2) * x[2], x[1]^2, x[2]^2, sqrt(2) * x[1] * x[2])
a <- c(1, 2); b <- c(-0.5, 3)
c(kernel = (1 + sum(a * b))^2, producto_en_espacio_ampliado = sum(phi(a) * phi(b)))

svm_pol <- svm(y ~ ., data = d_nl, kernel = "polynomial", degree = 3,
               gamma = 1, coef0 = 1, cost = 1)
svm_rad <- svm(y ~ ., data = d_nl, kernel = "radial", gamma = 0.5, cost = 1)

plot_svm(svm_pol, d_nl, "Kernel polinomial, d = 3") +
plot_svm(svm_rad, d_nl, "Kernel radial, γ = 0.5") + plot_layout(guides = "collect")

aj_rad <- svm(y ~ ., data = d_nl, kernel = "radial", gamma = 0.5, cost = 1, scale = FALSE)

K <- \(a, b) exp(-0.5 * sum((a - b)^2))                      # kernel radial con γ = 0.5
f_mano <- \(x) sum(aj_rad$coefs[, 1] * apply(aj_rad$SV, 1, K, b = x)) - aj_rad$rho

nuevos <- as.matrix(d_nl[1:6, c("x1", "x2")])
cbind(a_mano = apply(nuevos, 1, f_mano),
      e1071  = attr(predict(aj_rad, d_nl[1:6, ], decision.values = TRUE), "decision.values")[, 1])

heart <- read_csv("https://www.statlearning.com/s/Heart.csv", show_col_types = FALSE) |>
  select(-1) |>                              # primera columna = índice
  drop_na() |>                               # quita los 6 casos con NA (Ca, Thal)
  mutate(across(where(is.character), as.factor))
dim(heart)

set.seed(1)
idx      <- sample(nrow(heart), 207)         # 207 entrenamiento / 90 prueba, como en el libro
heart_tr <- heart[idx, ]
heart_te <- heart[-idx, ]

lda_fit <- MASS::lda(AHD ~ ., data = heart_tr)
svc     <- svm(AHD ~ ., data = heart_tr, kernel = "linear", cost = 1)
svm_g   <- map(c(`γ = 1e-3` = 1e-3, `γ = 1e-2` = 1e-2, `γ = 1e-1` = 1e-1),
               \(g) svm(AHD ~ ., data = heart_tr, kernel = "radial", gamma = g, cost = 1))

roc_de <- \(score, datos) roc(datos$AHD, score, levels = c("No", "Yes"),
                              direction = "<", quiet = TRUE)

rocs <- \(datos) list(
  izq = list(SVC = roc_de(score_svm(svc, datos), datos),
             LDA = roc_de(predict(lda_fit, datos)$posterior[, "Yes"], datos)),
  der = c(list(SVC = roc_de(score_svm(svc, datos), datos)),
          map(svm_g, \(m) roc_de(score_svm(m, datos), datos))))

tr <- rocs(heart_tr); te <- rocs(heart_te)
(plot_roc(tr$izq, "Entrenamiento: SVC vs. LDA") + plot_roc(tr$der, "Entrenamiento: kernel radial")) /
(plot_roc(te$izq, "Prueba: SVC vs. LDA")        + plot_roc(te$der, "Prueba: kernel radial"))

aj_svc <- svm(y ~ ., data = d_sol, kernel = "linear", cost = 1, scale = FALSE)
aj_glm <- glm(y ~ x1 + x2, data = d_sol, family = binomial)

cs <- coef_svm(aj_svc); cg <- coef(aj_glm)
fronteras <- tibble(modelo     = c("SVC (cost = 1)", "Logística"),
                    intercepto = c(-cs$beta0 / cs$beta[2], -cg[[1]] / cg[[3]]),
                    pendiente  = c(-cs$beta[1] / cs$beta[2], -cg[[2]] / cg[[3]]))

ggplot() +
  geom_point(data = d_sol, aes(x1, x2, colour = y), size = 2.5) +
  geom_abline(data = fronteras, aes(intercept = intercepto, slope = pendiente, linetype = modelo),
              linewidth = 1) +
  scale_colour_manual(values = pal) + coord_fixed()

set.seed(1)
d_r <- tibble(x = sort(runif(100, 0, 10)), y = sin(x) + rnorm(100, 0, 0.3))
svr <- svm(y ~ x, data = d_r, kernel = "radial", cost = 10, epsilon = 0.2)

d_r |> mutate(ajuste = predict(svr)) |>
  ggplot(aes(x, y)) + geom_point(alpha = .6) +
  geom_line(aes(y = ajuste), colour = pal[2], linewidth = 1.2) +
  geom_line(aes(y = ajuste + 0.2), colour = pal[2], linetype = "dashed") +
  geom_line(aes(y = ajuste - 0.2), colour = pal[2], linetype = "dashed")

set.seed(1)
x <- matrix(rnorm(20 * 2), ncol = 2)
y <- c(rep(-1, 10), rep(1, 10))
x[y == 1, ] <- x[y == 1, ] + 1
dat <- tibble(x1 = x[, 1], x2 = x[, 2], y = factor(y))

ggplot(dat, aes(x1, x2, colour = y)) + geom_point(size = 3) +
  scale_colour_manual(values = pal) + coord_fixed()

svmfit <- svm(y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)
plot_svm(svmfit, dat, "cost = 10")

svmfit$index          # índices de los vectores de soporte

summary(svmfit)

svmfit01 <- svm(y ~ ., data = dat, kernel = "linear", cost = 0.1, scale = FALSE)
plot_svm(svmfit01, dat, "cost = 0.1")

svmfit01$index

coef_svm(svmfit)

coef_svm(svmfit01)

set.seed(1)
tune_out <- tune(svm, y ~ ., data = dat, kernel = "linear",
                 ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))
summary(tune_out)

tune_out$performances |>
  ggplot(aes(cost, error)) +
  geom_line() + geom_point(size = 3) +
  geom_errorbar(aes(ymin = pmax(0, error - dispersion), ymax = error + dispersion), width = .1, alpha = .4) +
  scale_x_log10() +
  labs(x = "cost (escala log)", y = "Error de validación cruzada (10 pliegues)")

mejor <- tune_out$best.model

xtest <- matrix(rnorm(20 * 2), ncol = 2)
ytest <- sample(c(-1, 1), 20, replace = TRUE)
xtest[ytest == 1, ] <- xtest[ytest == 1, ] + 1
testdat <- tibble(x1 = xtest[, 1], x2 = xtest[, 2], y = factor(ytest))

ypred <- predict(mejor, testdat)
table(predicho = ypred, real = testdat$y)

mean(ypred == testdat$y)          # exactitud

# ¿y con cost = 0.01?
ypred_001 <- predict(svm(y ~ ., data = dat, kernel = "linear", cost = 0.01, scale = FALSE), testdat)
table(predicho = ypred_001, real = testdat$y)

x[y == 1, ] <- x[y == 1, ] + 0.5
dat_s <- tibble(x1 = x[, 1], x2 = x[, 2], y = factor(y))

fit_1e5 <- svm(y ~ ., data = dat_s, kernel = "linear", cost = 1e5, scale = FALSE)
fit_1   <- svm(y ~ ., data = dat_s, kernel = "linear", cost = 1,   scale = FALSE)

plot_svm(fit_1e5, dat_s, sprintf("cost = 1e5: %d SV, margen = %.2f", nrow(fit_1e5$SV), coef_svm(fit_1e5)$margen)) +
plot_svm(fit_1,   dat_s, sprintf("cost = 1: %d SV, margen = %.2f",   nrow(fit_1$SV),   coef_svm(fit_1)$margen)) +
  plot_layout(guides = "collect")

set.seed(1)
x <- matrix(rnorm(200 * 2), ncol = 2)
x[1:100, ]   <- x[1:100, ] + 2
x[101:150, ] <- x[101:150, ] - 2
y <- c(rep(1, 150), rep(2, 50))
dat <- tibble(x1 = x[, 1], x2 = x[, 2], y = factor(y))

train <- sample(200, 100)
ggplot(dat, aes(x1, x2, colour = y)) + geom_point(size = 2) + scale_colour_manual(values = pal)

svm_c1   <- svm(y ~ ., data = dat[train, ], kernel = "radial", gamma = 1, cost = 1)
svm_c1e5 <- svm(y ~ ., data = dat[train, ], kernel = "radial", gamma = 1, cost = 1e5)

plot_svm(svm_c1,   dat[train, ], "gamma = 1, cost = 1") +
plot_svm(svm_c1e5, dat[train, ], "gamma = 1, cost = 1e5") + plot_layout(guides = "collect")

summary(svm_c1)

set.seed(1)
tune_rad <- tune(svm, y ~ ., data = dat[train, ], kernel = "radial",
                 ranges = list(cost = c(0.1, 1, 10, 100, 1000), gamma = c(0.5, 1, 2, 3, 4)))
tune_rad$best.parameters

tune_rad$best.performance

tune_rad$performances |>
  ggplot(aes(factor(gamma), factor(cost), fill = error)) +
  geom_tile(colour = "white") + geom_text(aes(label = round(error, 2))) +
  scale_fill_gradient(low = "white", high = pal[1]) +
  labs(x = "gamma", y = "cost", fill = "Error CV")

pred_te <- predict(tune_rad$best.model, newdata = dat[-train, ])
table(real = dat[-train, ]$y, predicho = pred_te)

mean(pred_te != dat[-train, ]$y)       # tasa de error (el libro reporta ≈ 12 %)

plot_conf(pred_te, dat[-train, ]$y, "Prueba")

svm_opt  <- svm(y ~ ., data = dat[train, ], kernel = "radial", gamma = 2,  cost = 1)
svm_flex <- svm(y ~ ., data = dat[train, ], kernel = "radial", gamma = 50, cost = 1)
modelos  <- list("γ = 2" = svm_opt, "γ = 50" = svm_flex)

# score_svm(): positivo = segunda clase  =>  direction = "<" con levels = c("1", "2")
rocs_lab <- \(datos) map(modelos, \(m) roc(datos$y, score_svm(m, datos),
                                           levels = levels(datos$y), direction = "<", quiet = TRUE))

plot_roc(rocs_lab(dat[train, ]),  "Entrenamiento") +
plot_roc(rocs_lab(dat[-train, ]), "Prueba")

set.seed(1)
x <- rbind(x, matrix(rnorm(50 * 2), ncol = 2))
y <- c(y, rep(0, 50))
x[y == 0, 2] <- x[y == 0, 2] + 2
dat3 <- tibble(x1 = x[, 1], x2 = x[, 2], y = factor(y))

svm3 <- svm(y ~ ., data = dat3, kernel = "radial", cost = 10, gamma = 1)
plot_svm(svm3, dat3, "SVM radial con 3 clases (uno contra uno)")

# Uno contra uno: e1071 ajusta choose(3, 2) = 3 clasificadores
attr(predict(svm3, dat3[1:4, ], decision.values = TRUE), "decision.values")

# Uno contra todos: K clasificadores binarios, gana el mayor puntaje
ova <- map(levels(dat3$y), \(k) {
  d <- mutate(dat3, y = factor(if_else(y == k, "si", "no"), levels = c("no", "si")))
  svm(y ~ ., data = d, kernel = "radial", cost = 10, gamma = 1)
})
puntajes <- map(ova, \(m) score_svm(m, dat3)) |> set_names(levels(dat3$y)) |> as_tibble()
pred_ova <- levels(dat3$y)[max.col(puntajes)]

c(exactitud_OvO = mean(predict(svm3, dat3) == dat3$y),
  exactitud_OvA = mean(pred_ova == dat3$y))

names(Khan)

dim(Khan$xtrain); dim(Khan$xtest)

table(entrenamiento = Khan$ytrain)

table(prueba = Khan$ytest)

out <- svm(Khan$xtrain, factor(Khan$ytrain), kernel = "linear", cost = 10)
summary(out)

table(ajustado = out$fitted, real = Khan$ytrain)      # entrenamiento: sin errores

pred_te <- predict(out, Khan$xtest)
table(predicho = pred_te, real = Khan$ytest)           # prueba

plot_conf(pred_te, factor(Khan$ytest), "Khan: prueba (cost = 10)")

set.seed(1)
tune_khan <- tune.svm(Khan$xtrain, factor(Khan$ytrain), kernel = "linear",
                      cost = c(0.001, 0.01, 0.1, 1, 10, 100))
tune_khan$performances |>
  ggplot(aes(cost, error)) + geom_line() + geom_point(size = 3) + scale_x_log10() +
  labs(x = "cost (escala log)", y = "Error de CV")

table(predicho = predict(tune_khan$best.model, Khan$xtest), real = Khan$ytest)