ISLR - Lab 9.6: Support Vector Machines

9.6.1 Support Vector Classifier

La libreria e1071 contiene muchos de modelos estadisticos y uno de ellos es svm() que vamos a utilizar para realizar un modelo maquina de vectores de soporte sin olvidar el parametro kernel=“linear”, también tenemos el parámetro cost que cuando es bajo nos permite ajustar mas vectores y cuando es muy alto permine menos vectores.

Generemos los datos para demostrar el uso de cost y comprobemos si son linealmente separables.

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
plot(x, col = (3 - y))

Vemos que no lo son. Ahora ajustaremos el modelo, notemos que para hacer la clasificación debemos codificar como factor la respuesta.

dat = data.frame(x = x, y = as.factor (y))
library (e1071)
svmfit = svm(
  y~.,
  data = dat ,
  kernel ="linear",
  cost = 10,
  scale = FALSE
)

el parmetro scale=FALSE indica que no se se haga la escala de cada variable hasta llegar a cero o la primer desviacion estandar

plot(svmfit , dat)

Lo que se encuentra en celeste sera asignado a -1 y lo lila sera asignado a 1, la frontera de decision entre las 2 clases es lineal por el parametro kernel que utilizamos. Los vectores son las x y podemos verlos con

svmfit$index

[1]  1  2  5  7 14 16 17

podemos obtener mas informacion sobre el modelo asi

summary (svmfit )


Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  10 

Number of Support Vectors:  7

 ( 4 3 )


Number of Classes:  2 

Levels: 
 -1 1

Ahora veamos que sucede al utilizar un cost mas alto

svmfit = svm(
  y~.,
  data = dat ,
  kernel ="linear",
  cost = 0.1,
  scale = FALSE
)
plot(svmfit , dat)

svmfit$index

 [1]  1  2  3  4  5  7  9 10 12 13 14 15 16 17 18 20

esta libreria tambien cuenta con la funcion tune() para realizar k-fold CV que por defecto toma k=10. A continuacion de como cofigurarlo utilizando distintos valores en cost.

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)))

y con summary() podemos consultar los errores de validacion cruzada

summary (tune.out)


Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:


- best performance: 0.1 

- Detailed performance results:

NA

vemos que cost = 0.1 es el mejor resultado y tune() por defecto lo almacena en best.model

bestmod =tune.out$best.model
summary (bestmod )


Call:
best.tune(method = svm, train.x = y ~ ., data = dat, ranges = list(cost = c(0.001, 
    0.01, 0.1, 1, 5, 10, 100)), kernel = "linear")


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  0.1 

Number of Support Vectors:  16

 ( 8 8 )


Number of Classes:  2 

Levels: 
 -1 1

luego podemoos utilizar la funcion predict() para clasificar un set de datos de prueba y utilizaremos el mejor modelo segun la validacion cruzada

xtest=matrix (rnorm (20*2) , ncol =2)
ytest=sample (c(-1,1) , 20, rep=TRUE)
xtest[ytest ==1 ,]= xtest[ytest ==1,] + 1
testdat =data.frame (x=xtest , y=as.factor (ytest))

ypred=predict (bestmod ,testdat )
table(predict =ypred , truth= testdat$y )

       truth
predict -1  1
     -1 11  1
     1   0  8

con esto conseguimos 19 clasificaciones correctas pero que sucede si utilizamos cost=0.01

svmfit = svm(
  y~.,
  data = dat ,
  kernel ="linear",
  cost = .01,
  scale = FALSE
)
ypred = predict (svmfit , testdat)
table(predict = ypred , truth = testdat$y)

       truth
predict -1  1
     -1 11  2
     1   0  7

y vemos que tenemos un resultado correcto menos.

Ahora asumamos que las 2 clases son linealmente separables, configuremos nuestro set de datos y ajustemos el modelo con un valor cost muy alto

dat = data.frame(x = x, y = as.factor (y))
svmfit = svm(y~.,
             data = dat ,
             kernel ="linear",
             cost = 1e5)
summary (svmfit)


Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1e+05)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  1e+05 

Number of Support Vectors:  7

 ( 4 3 )


Number of Classes:  2 

Levels: 
 -1 1

plot(svmfit , dat)

solo obtenemos 3 vectores de soporte porque en la figura vemos que las observaciones que no son vectores de soporte estan muy cerca de la frontera de decision, por ello asumimos que tendremos malos resultados si probamos un dataset de prueba.

Ahora configuremos un cost mas bajo

svmfit = svm(y~.,
             data = dat ,
             kernel ="linear",
             cost = 1)
summary (svmfit)


Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  1 

Number of Support Vectors:  11

 ( 6 5 )


Number of Classes:  2 

Levels: 
 -1 1

plot(svmfit , dat)

con cost=1 solo fallamos en una observacion y obtuvimos 7 vectores. Parece que este modelo tendra mejores resultados al hacer una comprobacion.

9.6.2 Support Vector Machine

Para un maquina de soporte de verctores no lineal utilizamos 2 nuevos parametros kernel=“polynomimal” y kernel=“radial”, tambien tenemos degree para indicar los grados de la forma polinomica y gamma para los kernel radiales.

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 = data.frame(x = x, y = as.factor (y))

al graficar vemos claramente que no es lineal la frontera

plot(x, col=y)

dividimos en train y test y ajustamos un melo radial con gamma=1

train = sample (200 , 100)
svmfit = svm(
  y~.,
  data = dat [train , ],
  kernel ="radial",
  gamma = 1,
  cost = 1
)
plot(svmfit , dat[train , ])


summary (svmfit)


Call:
svm(formula = y ~ ., data = dat[train, ], kernel = "radial", gamma = 1, cost = 1)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 

Number of Support Vectors:  37

 ( 17 20 )


Number of Classes:  2 

Levels: 
 1 2

en la grafica podemos ver que tenemos muchos errores de clasificacion, si aumentamos el costo podemos reducir errores pero esto viene a costas de una frontera mas irregular con la que corremos el riesgo de tener sobreajuste.

svmfit = svm(
  y~.,
  data = dat [train , ],
  kernel ="radial",
  gamma = 1,
  cost = 1e5
)
plot(svmfit , dat [train , ])

esto tambien lo podemos combinar con tune() para seleccionar la mejor combinacion de cost y gamma

set.seed (1)
tune.out = 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)
  )
)
summary (tune.out)


Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:


- best performance: 0.12 

- Detailed performance results:

NA

los mejores resultados fueron con cost=1 y gamma=2, probemos con el set de prueba

table(true = dat[-train ,"y"],
      pred = predict (tune.out$best.model , newdata = dat[-train , ]))

tenemos un 10% de error.

9.6.3 ROC Curves

Utilizamos el paquete ROCR para para generar curvas ROC, primero creamos una funcion para graficar las curvas dado un vector con las observaciones. En una maquin a de soporte de vectores la clasificacion de la etiqueta se hace segun el signo del factor.

library (ROCR)

Loading required package: gplots

Attaching package: 㤼㸱gplots㤼㸲

The following object is masked from 㤼㸱package:stats㤼㸲:

    lowess

rocplot = function (pred , truth , ...) {
  predob = prediction (pred , truth)
  perf = performance (predob ,"tpr","fpr")
  plot(perf , ...)
}

svmfit.opt = svm(
  y~.,
  data = dat[train , ],
  kernel ="radial",
  gamma = 2,
  cost = 1,
  decision.values = T
)

fitted =attributes (predict (svmfit.opt ,dat[train ,], decision.values =TRUE))$decision.values

ahora graficamos el ROC

# par(mfrow =c(1,2))
# rocplot (fitted ,dat[train,"y"], main="Training Data")

parece estar haciendo un prediccion decente pero podemos aumentar gamma para mejorar los resultados

svmfit.flex = svm (
  y~.,
  data = dat[train , ],
  kernel ="radial",
  gamma = 50,
  cost = 1,
  decision.values = T
)
fitted = attributes (predict (svmfit.flex , dat[train , ], decision.values =
                                T))$decision.values 
# rocplot(fitted , dat [train ,"y"], add = T, col ="red")

luego probamo el modelo con set de prueba el modelo con gamma=2 parece dar los mejores resultados

fitted = attributes (predict (svmfit.opt , dat[-train , ], decision.values =
                                T))$decision.values
# rocplot (fitted , dat [-train ,"y"], main ="Test Data")
fitted = attributes (predict (svmfit.flex , dat[-train , ], decision.values =
                                T))$decision.values

str(fitted)

 num [1:100, 1] -0.946 -0.505 -0.968 -0.726 -0.589 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:100] "12" "14" "15" "16" ...
  ..$ : chr "2/1"

# rocplot (fitted , dat [-train ,"y"], add = T, col ="red")

SVM with Multiple Classes

Si la respuesta que buscamos es un factor con mas de 2 valores svm() realizará una clasificación multiclase

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
dat = data.frame(x = x, y = as.factor (y))
par(mfrow = c(1, 1))
plot(x, col = (y + 1))


svmfit = svm(
  y~.,
  data = dat ,
  kernel ="radial",
  cost = 10,
  gamma = 1
)
plot(svmfit , dat)

la funcion svm() tambien puede hacer maquinas de soporte de vectores de regresion si la variable es un numerica y no un factor.

9.6.5 Application to Gene Expression Data

usaremos el set de datos Khan acerca de 4 tejidos cancerigenos.

library (ISLR)
names(Khan)

[1] "xtrain" "xtest"  "ytrain" "ytest"

dim(Khan$xtrain)

[1]   63 2308

dim(Khan$xtest)

[1]   20 2308

dim(Khan$ytrain)

NULL

dim(Khan$ytest)

NULL

table(Khan$ytrain )


 1  2  3  4 
 8 23 12 20

table(Khan$ytest )


1 2 3 4 
3 6 6 5

usaremos un modelo lineal por la gran cantidad de observaciones que tenemos


dat=data.frame(x=Khan$xtrain , y=as.factor( Khan$ytrain ))
out=svm(y~., data=dat , kernel ="linear",cost =10)
summary (out)


Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  10 

Number of Support Vectors:  58

 ( 20 20 11 7 )


Number of Classes:  4 

Levels: 
 1 2 3 4

table(out$fitted , dat$y)

   
     1  2  3  4
  1  8  0  0  0
  2  0 23  0  0
  3  0  0 12  0
  4  0  0  0 20

###libera memoria

gc()

          used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells 2012002 107.5    3422272 182.8  3422272 182.8
Vcells 5298944  40.5   18104948 138.2 19805804 151.2

LS0tDQp0aXRsZTogJ0lTTFIgLSBMYWIgOS42OiBTdXBwb3J0IFZlY3RvciBNYWNoaW5lcycNCm91dHB1dDoNCiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdA0KICBodG1sX2RvY3VtZW50Og0KICAgIGRmX3ByaW50OiBwYWdlZA0KICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQNCi0tLQ0KDQojIyMgOS42LjEgU3VwcG9ydCBWZWN0b3IgQ2xhc3NpZmllcg0KDQpMYSBsaWJyZXJpYSAqZTEwNzEqIGNvbnRpZW5lIG11Y2hvcyBkZSBtb2RlbG9zIGVzdGFkaXN0aWNvcyB5IHVubyBkZSBlbGxvcyBlcyAqc3ZtKCkqIHF1ZSB2YW1vcyBhIHV0aWxpemFyIHBhcmEgcmVhbGl6YXIgdW4gbW9kZWxvICptYXF1aW5hIGRlIHZlY3RvcmVzIGRlIHNvcG9ydGUqIHNpbiBvbHZpZGFyIGVsIHBhcmFtZXRybyAqa2VybmVsPSJsaW5lYXIiKiwgdGFtYmnDqW4gdGVuZW1vcyBlbCBwYXLDoW1ldHJvICpjb3N0KiBxdWUgY3VhbmRvIGVzIGJham8gbm9zIHBlcm1pdGUgYWp1c3RhciBtYXMgdmVjdG9yZXMgeSBjdWFuZG8gZXMgbXV5IGFsdG8gcGVybWluZSBtZW5vcyB2ZWN0b3Jlcy4NCg0KR2VuZXJlbW9zIGxvcyBkYXRvcyBwYXJhIGRlbW9zdHJhciBlbCB1c28gZGUgKmNvc3QqIHkgY29tcHJvYmVtb3Mgc2kgc29uIGxpbmVhbG1lbnRlIHNlcGFyYWJsZXMuDQoNCmBgYHtyfQ0Kc2V0LnNlZWQgKDEpDQp4ID0gbWF0cml4IChybm9ybSAoMjAgKiAyKSAsIG5jb2wgPSAyKQ0KeSA9IGMocmVwICgtMSwgMTApICwgcmVwICgxICwgMTApKQ0KeFt5ID09IDEgLCBdID0geFt5ID09IDEsIF0gKyAxDQpwbG90KHgsIGNvbCA9ICgzIC0geSkpDQoNCmBgYA0KDQpWZW1vcyBxdWUgbm8gbG8gc29uLiBBaG9yYSBhanVzdGFyZW1vcyBlbCBtb2RlbG8sIG5vdGVtb3MgcXVlIHBhcmEgaGFjZXIgbGEgY2xhc2lmaWNhY2nDs24gZGViZW1vcyBjb2RpZmljYXIgY29tbyBmYWN0b3IgbGEgcmVzcHVlc3RhLg0KDQpgYGB7cn0NCmRhdCA9IGRhdGEuZnJhbWUoeCA9IHgsIHkgPSBhcy5mYWN0b3IgKHkpKQ0KbGlicmFyeSAoZTEwNzEpDQpzdm1maXQgPSBzdm0oDQogIHl+LiwNCiAgZGF0YSA9IGRhdCAsDQogIGtlcm5lbCA9ImxpbmVhciIsDQogIGNvc3QgPSAxMCwNCiAgc2NhbGUgPSBGQUxTRQ0KKQ0KYGBgDQoNCmVsIHBhcm1ldHJvICpzY2FsZT1GQUxTRSogaW5kaWNhIHF1ZSBubyBzZSBzZSBoYWdhIGxhIGVzY2FsYSBkZSBjYWRhIHZhcmlhYmxlIGhhc3RhIGxsZWdhciBhIGNlcm8gbyBsYSBwcmltZXIgZGVzdmlhY2lvbiBlc3RhbmRhcg0KDQpgYGB7cn0NCnBsb3Qoc3ZtZml0ICwgZGF0KQ0KYGBgDQoNCkxvIHF1ZSBzZSBlbmN1ZW50cmEgZW4gY2VsZXN0ZSBzZXJhIGFzaWduYWRvIGEgLTEgeSBsbyBsaWxhIHNlcmEgYXNpZ25hZG8gYSAxLCBsYSBmcm9udGVyYSBkZSBkZWNpc2lvbiBlbnRyZSBsYXMgMiBjbGFzZXMgZXMgbGluZWFsIHBvciBlbCBwYXJhbWV0cm8gKmtlcm5lbCogcXVlIHV0aWxpemFtb3MuIExvcyB2ZWN0b3JlcyBzb24gbGFzIHggeSBwb2RlbW9zIHZlcmxvcyBjb24gDQoNCmBgYHtyfQ0Kc3ZtZml0JGluZGV4DQpgYGANCg0KcG9kZW1vcyBvYnRlbmVyIG1hcyBpbmZvcm1hY2lvbiBzb2JyZSBlbCBtb2RlbG8gYXNpDQoNCmBgYHtyfQ0Kc3VtbWFyeSAoc3ZtZml0ICkNCmBgYA0KDQpBaG9yYSB2ZWFtb3MgcXVlIHN1Y2VkZSBhbCB1dGlsaXphciB1biAqY29zdCogbWFzIGFsdG8NCg0KYGBge3J9DQpzdm1maXQgPSBzdm0oDQogIHl+LiwNCiAgZGF0YSA9IGRhdCAsDQogIGtlcm5lbCA9ImxpbmVhciIsDQogIGNvc3QgPSAwLjEsDQogIHNjYWxlID0gRkFMU0UNCikNCnBsb3Qoc3ZtZml0ICwgZGF0KQ0Kc3ZtZml0JGluZGV4DQpgYGANCg0KZXN0YSBsaWJyZXJpYSB0YW1iaWVuIGN1ZW50YSBjb24gbGEgZnVuY2lvbiAqdHVuZSgpKiBwYXJhIHJlYWxpemFyICprLWZvbGQgQ1YqIHF1ZSBwb3IgZGVmZWN0byB0b21hICprPTEwKi4gQSBjb250aW51YWNpb24gZGUgY29tbyBjb2ZpZ3VyYXJsbyB1dGlsaXphbmRvIGRpc3RpbnRvcyB2YWxvcmVzIGVuICpjb3N0Ki4NCg0KYGBge3J9DQpzZXQuc2VlZCAoMSkNCnR1bmUub3V0ID0gdHVuZShzdm0gLA0KICAgICAgICAgICAgICAgIHl+LiwNCiAgICAgICAgICAgICAgICBkYXRhID0gZGF0ICwNCiAgICAgICAgICAgICAgICBrZXJuZWwgPSJsaW5lYXIiLA0KICAgICAgICAgICAgICAgIHJhbmdlcyA9IGxpc3QoY29zdCA9IGMoMC4wMDEgLCAwLjAxLCAwLjEsIDEsIDUsIDEwLCAxMDApKSkNCmBgYA0KDQp5IGNvbiAqc3VtbWFyeSgpKiBwb2RlbW9zIGNvbnN1bHRhciBsb3MgZXJyb3JlcyBkZSB2YWxpZGFjaW9uIGNydXphZGENCg0KYGBge3J9DQpzdW1tYXJ5ICh0dW5lLm91dCkNCmBgYA0KDQp2ZW1vcyBxdWUgY29zdCA9IDAuMSAgZXMgZWwgbWVqb3IgcmVzdWx0YWRvIHkgKnR1bmUoKSogcG9yIGRlZmVjdG8gbG8gYWxtYWNlbmEgZW4gKmJlc3QubW9kZWwqDQoNCmBgYHtyfQ0KYmVzdG1vZCA9dHVuZS5vdXQkYmVzdC5tb2RlbA0Kc3VtbWFyeSAoYmVzdG1vZCApDQpgYGANCg0KbHVlZ28gcG9kZW1vb3MgdXRpbGl6YXIgbGEgZnVuY2lvbiAqcHJlZGljdCgpKiBwYXJhIGNsYXNpZmljYXIgdW4gc2V0IGRlIGRhdG9zIGRlIHBydWViYSB5IHV0aWxpemFyZW1vcyBlbCBtZWpvciBtb2RlbG8gc2VndW4gbGEgdmFsaWRhY2lvbiBjcnV6YWRhDQoNCmBgYHtyfQ0KeHRlc3Q9bWF0cml4IChybm9ybSAoMjAqMikgLCBuY29sID0yKQ0KeXRlc3Q9c2FtcGxlIChjKC0xLDEpICwgMjAsIHJlcD1UUlVFKQ0KeHRlc3RbeXRlc3QgPT0xICxdPSB4dGVzdFt5dGVzdCA9PTEsXSArIDENCnRlc3RkYXQgPWRhdGEuZnJhbWUgKHg9eHRlc3QgLCB5PWFzLmZhY3RvciAoeXRlc3QpKQ0KDQp5cHJlZD1wcmVkaWN0IChiZXN0bW9kICx0ZXN0ZGF0ICkNCnRhYmxlKHByZWRpY3QgPXlwcmVkICwgdHJ1dGg9IHRlc3RkYXQkeSApDQpgYGANCg0KY29uIGVzdG8gY29uc2VndWltb3MgMTkgY2xhc2lmaWNhY2lvbmVzIGNvcnJlY3RhcyBwZXJvIHF1ZSBzdWNlZGUgc2kgdXRpbGl6YW1vcyAqY29zdD0wLjAxKg0KDQpgYGB7cn0NCnN2bWZpdCA9IHN2bSgNCiAgeX4uLA0KICBkYXRhID0gZGF0ICwNCiAga2VybmVsID0ibGluZWFyIiwNCiAgY29zdCA9IC4wMSwNCiAgc2NhbGUgPSBGQUxTRQ0KKQ0KeXByZWQgPSBwcmVkaWN0IChzdm1maXQgLCB0ZXN0ZGF0KQ0KdGFibGUocHJlZGljdCA9IHlwcmVkICwgdHJ1dGggPSB0ZXN0ZGF0JHkpDQpgYGANCg0KeSB2ZW1vcyBxdWUgdGVuZW1vcyB1biByZXN1bHRhZG8gY29ycmVjdG8gbWVub3MuDQoNCkFob3JhIGFzdW1hbW9zIHF1ZSBsYXMgMiBjbGFzZXMgc29uIGxpbmVhbG1lbnRlIHNlcGFyYWJsZXMsIGNvbmZpZ3VyZW1vcyBudWVzdHJvIHNldCBkZSBkYXRvcyB5IGFqdXN0ZW1vcyBlbCBtb2RlbG8gY29uIHVuIHZhbG9yICpjb3N0KiBtdXkgYWx0bw0KDQpgYGB7cn0NCmRhdCA9IGRhdGEuZnJhbWUoeCA9IHgsIHkgPSBhcy5mYWN0b3IgKHkpKQ0Kc3ZtZml0ID0gc3ZtKHl+LiwNCiAgICAgICAgICAgICBkYXRhID0gZGF0ICwNCiAgICAgICAgICAgICBrZXJuZWwgPSJsaW5lYXIiLA0KICAgICAgICAgICAgIGNvc3QgPSAxZTUpDQpzdW1tYXJ5IChzdm1maXQpDQpwbG90KHN2bWZpdCAsIGRhdCkNCmBgYA0KDQpzb2xvIG9idGVuZW1vcyAzIHZlY3RvcmVzIGRlIHNvcG9ydGUgcG9ycXVlIGVuIGxhIGZpZ3VyYSB2ZW1vcyBxdWUgbGFzIG9ic2VydmFjaW9uZXMgcXVlIG5vIHNvbiB2ZWN0b3JlcyBkZSBzb3BvcnRlIGVzdGFuIG11eSBjZXJjYSBkZSBsYSBmcm9udGVyYSBkZSBkZWNpc2lvbiwgcG9yIGVsbG8gYXN1bWltb3MgcXVlIHRlbmRyZW1vcyBtYWxvcyByZXN1bHRhZG9zIHNpIHByb2JhbW9zIHVuIGRhdGFzZXQgZGUgcHJ1ZWJhLg0KDQpBaG9yYSBjb25maWd1cmVtb3MgdW4gKmNvc3QqIG1hcyBiYWpvDQoNCmBgYHtyfQ0Kc3ZtZml0ID0gc3ZtKHl+LiwNCiAgICAgICAgICAgICBkYXRhID0gZGF0ICwNCiAgICAgICAgICAgICBrZXJuZWwgPSJsaW5lYXIiLA0KICAgICAgICAgICAgIGNvc3QgPSAxKQ0Kc3VtbWFyeSAoc3ZtZml0KQ0KcGxvdChzdm1maXQgLCBkYXQpDQpgYGANCg0KY29uICpjb3N0PTEqIHNvbG8gZmFsbGFtb3MgZW4gdW5hIG9ic2VydmFjaW9uIHkgb2J0dXZpbW9zIDcgdmVjdG9yZXMuIFBhcmVjZSBxdWUgZXN0ZSBtb2RlbG8gdGVuZHJhIG1lam9yZXMgcmVzdWx0YWRvcyBhbCBoYWNlciB1bmEgY29tcHJvYmFjaW9uLg0KDQojIyMgOS42LjIgU3VwcG9ydCBWZWN0b3IgTWFjaGluZQ0KDQpQYXJhIHVuIG1hcXVpbmEgZGUgc29wb3J0ZSBkZSB2ZXJjdG9yZXMgbm8gbGluZWFsIHV0aWxpemFtb3MgMiBudWV2b3MgcGFyYW1ldHJvcyAqa2VybmVsPSJwb2x5bm9taW1hbCIqIHkgKmtlcm5lbD0icmFkaWFsIiosIHRhbWJpZW4gdGVuZW1vcyAqZGVncmVlKiBwYXJhIGluZGljYXIgbG9zIGdyYWRvcyBkZSBsYSBmb3JtYSBwb2xpbm9taWNhIHkgKmdhbW1hKiBwYXJhIGxvcyBrZXJuZWwgcmFkaWFsZXMuDQoNCmBgYHtyfQ0Kc2V0LnNlZWQgKDEpDQp4ID0gbWF0cml4IChybm9ybSAoMjAwICogMikgLCBuY29sID0gMikNCnhbMToxMDAgLCBdID0geFsxOjEwMCAsIF0gKyAyDQp4WzEwMToxNTAgLCBdID0geFsxMDE6MTUwICwgXSAtIDINCnkgPSBjKHJlcCAoMSAsIDE1MCkgLCByZXAgKDIgLCA1MCkpDQpkYXQgPSBkYXRhLmZyYW1lKHggPSB4LCB5ID0gYXMuZmFjdG9yICh5KSkNCmBgYA0KDQphbCBncmFmaWNhciB2ZW1vcyBjbGFyYW1lbnRlIHF1ZSBubyBlcyBsaW5lYWwgbGEgZnJvbnRlcmENCg0KYGBge3J9DQpwbG90KHgsIGNvbD15KQ0KYGBgDQoNCmRpdmlkaW1vcyBlbiAqdHJhaW4qIHkgKnRlc3QqIHkgYWp1c3RhbW9zIHVuIG1lbG8gKnJhZGlhbCogY29uICpnYW1tYT0xKiANCg0KYGBge3J9DQp0cmFpbiA9IHNhbXBsZSAoMjAwICwgMTAwKQ0Kc3ZtZml0ID0gc3ZtKA0KICB5fi4sDQogIGRhdGEgPSBkYXQgW3RyYWluICwgXSwNCiAga2VybmVsID0icmFkaWFsIiwNCiAgZ2FtbWEgPSAxLA0KICBjb3N0ID0gMQ0KKQ0KcGxvdChzdm1maXQgLCBkYXRbdHJhaW4gLCBdKQ0KDQpzdW1tYXJ5IChzdm1maXQpDQpgYGANCg0KZW4gbGEgZ3JhZmljYSBwb2RlbW9zIHZlciBxdWUgdGVuZW1vcyBtdWNob3MgZXJyb3JlcyBkZSBjbGFzaWZpY2FjaW9uLCBzaSBhdW1lbnRhbW9zIGVsICpjb3N0byogcG9kZW1vcyByZWR1Y2lyIGVycm9yZXMgcGVybyBlc3RvIHZpZW5lIGEgY29zdGFzIGRlIHVuYSBmcm9udGVyYSBtYXMgaXJyZWd1bGFyIGNvbiBsYSBxdWUgY29ycmVtb3MgZWwgcmllc2dvIGRlIHRlbmVyIHNvYnJlYWp1c3RlLg0KDQpgYGB7cn0NCnN2bWZpdCA9IHN2bSgNCiAgeX4uLA0KICBkYXRhID0gZGF0IFt0cmFpbiAsIF0sDQogIGtlcm5lbCA9InJhZGlhbCIsDQogIGdhbW1hID0gMSwNCiAgY29zdCA9IDFlNQ0KKQ0KcGxvdChzdm1maXQgLCBkYXQgW3RyYWluICwgXSkNCmBgYA0KDQplc3RvIHRhbWJpZW4gbG8gcG9kZW1vcyBjb21iaW5hciBjb24gKnR1bmUoKSogcGFyYSBzZWxlY2Npb25hciBsYSBtZWpvciBjb21iaW5hY2lvbiBkZSAqY29zdCogeSAqZ2FtbWEqDQoNCmBgYHtyfQ0Kc2V0LnNlZWQgKDEpDQp0dW5lLm91dCA9IHR1bmUoDQogIHN2bSAsDQogIHl+LiwNCiAgZGF0YSA9IGRhdFt0cmFpbiAsIF0sDQogIGtlcm5lbCA9InJhZGlhbCIsDQogIHJhbmdlcyA9IGxpc3QoDQogICAgY29zdCA9IGMoMC4xICwgMSAsIDEwICwgMTAwICwgMTAwMCksDQogICAgZ2FtbWEgPSBjKDAuNSwgMSwgMiwgMywgNCkNCiAgKQ0KKQ0Kc3VtbWFyeSAodHVuZS5vdXQpDQpgYGANCg0KbG9zIG1lam9yZXMgcmVzdWx0YWRvcyBmdWVyb24gY29uICpjb3N0PTEqIHkgKmdhbW1hPTIqLCBwcm9iZW1vcyBjb24gZWwgc2V0IGRlIHBydWViYQ0KDQpgYGB7cn0NCnRhYmxlKHRydWUgPSBkYXRbLXRyYWluICwieSJdLA0KICAgICAgcHJlZCA9IHByZWRpY3QgKHR1bmUub3V0JGJlc3QubW9kZWwgLCBuZXdkYXRhID0gZGF0Wy10cmFpbiAsIF0pKQ0KYGBgDQoNCnRlbmVtb3MgdW4gMTAlIGRlIGVycm9yLg0KDQojIyMgOS42LjMgUk9DIEN1cnZlcw0KDQpVdGlsaXphbW9zIGVsIHBhcXVldGUgUk9DUiBwYXJhIHBhcmEgZ2VuZXJhciBjdXJ2YXMgUk9DLCBwcmltZXJvIGNyZWFtb3MgdW5hIGZ1bmNpb24gcGFyYSBncmFmaWNhciBsYXMgY3VydmFzIGRhZG8gdW4gdmVjdG9yIGNvbiBsYXMgb2JzZXJ2YWNpb25lcy4gIEVuIHVuYSBtYXF1aW4gYSBkZSBzb3BvcnRlIGRlIHZlY3RvcmVzIGxhIGNsYXNpZmljYWNpb24gZGUgbGEgZXRpcXVldGEgc2UgaGFjZSBzZWd1biBlbCBzaWdubyBkZWwgZmFjdG9yLg0KDQpgYGB7cn0NCmxpYnJhcnkgKFJPQ1IpDQpyb2NwbG90ID0gZnVuY3Rpb24gKHByZWQgLCB0cnV0aCAsIC4uLikgew0KICBwcmVkb2IgPSBwcmVkaWN0aW9uIChwcmVkICwgdHJ1dGgpDQogIHBlcmYgPSBwZXJmb3JtYW5jZSAocHJlZG9iICwidHByIiwiZnByIikNCiAgcGxvdChwZXJmICwgLi4uKQ0KfQ0KDQpzdm1maXQub3B0ID0gc3ZtKA0KICB5fi4sDQogIGRhdGEgPSBkYXRbdHJhaW4gLCBdLA0KICBrZXJuZWwgPSJyYWRpYWwiLA0KICBnYW1tYSA9IDIsDQogIGNvc3QgPSAxLA0KICBkZWNpc2lvbi52YWx1ZXMgPSBUDQopDQoNCmZpdHRlZCA9YXR0cmlidXRlcyAocHJlZGljdCAoc3ZtZml0Lm9wdCAsZGF0W3RyYWluICxdLCBkZWNpc2lvbi52YWx1ZXMgPVRSVUUpKSRkZWNpc2lvbi52YWx1ZXMNCmBgYA0KDQphaG9yYSBncmFmaWNhbW9zIGVsIFJPQw0KDQpgYGB7cn0NCiMgcGFyKG1mcm93ID1jKDEsMikpDQojIHJvY3Bsb3QgKGZpdHRlZCAsZGF0W3RyYWluLCJ5Il0sIG1haW49IlRyYWluaW5nIERhdGEiKQ0KYGBgDQoNCnBhcmVjZSBlc3RhciBoYWNpZW5kbyB1biBwcmVkaWNjaW9uIGRlY2VudGUgcGVybyBwb2RlbW9zIGF1bWVudGFyICpnYW1tYSogcGFyYSBtZWpvcmFyIGxvcyByZXN1bHRhZG9zDQoNCmBgYHtyfQ0Kc3ZtZml0LmZsZXggPSBzdm0gKA0KICB5fi4sDQogIGRhdGEgPSBkYXRbdHJhaW4gLCBdLA0KICBrZXJuZWwgPSJyYWRpYWwiLA0KICBnYW1tYSA9IDUwLA0KICBjb3N0ID0gMSwNCiAgZGVjaXNpb24udmFsdWVzID0gVA0KKQ0KZml0dGVkID0gYXR0cmlidXRlcyAocHJlZGljdCAoc3ZtZml0LmZsZXggLCBkYXRbdHJhaW4gLCBdLCBkZWNpc2lvbi52YWx1ZXMgPQ0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBUKSkkZGVjaXNpb24udmFsdWVzIA0KIyByb2NwbG90KGZpdHRlZCAsIGRhdCBbdHJhaW4gLCJ5Il0sIGFkZCA9IFQsIGNvbCA9InJlZCIpDQoNCmBgYA0KDQpsdWVnbyBwcm9iYW1vIGVsIG1vZGVsbyBjb24gc2V0IGRlIHBydWViYSBlbCBtb2RlbG8gY29uICpnYW1tYT0yKiBwYXJlY2UgZGFyIGxvcyBtZWpvcmVzIHJlc3VsdGFkb3MNCg0KYGBge3J9DQpmaXR0ZWQgPSBhdHRyaWJ1dGVzIChwcmVkaWN0IChzdm1maXQub3B0ICwgZGF0Wy10cmFpbiAsIF0sIGRlY2lzaW9uLnZhbHVlcyA9DQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFQpKSRkZWNpc2lvbi52YWx1ZXMNCiMgcm9jcGxvdCAoZml0dGVkICwgZGF0IFstdHJhaW4gLCJ5Il0sIG1haW4gPSJUZXN0IERhdGEiKQ0KZml0dGVkID0gYXR0cmlidXRlcyAocHJlZGljdCAoc3ZtZml0LmZsZXggLCBkYXRbLXRyYWluICwgXSwgZGVjaXNpb24udmFsdWVzID0NCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgVCkpJGRlY2lzaW9uLnZhbHVlcw0KDQpzdHIoZml0dGVkKQ0KIyByb2NwbG90IChmaXR0ZWQgLCBkYXQgWy10cmFpbiAsInkiXSwgYWRkID0gVCwgY29sID0icmVkIikNCmBgYA0KDQojIyMgU1ZNIHdpdGggTXVsdGlwbGUgQ2xhc3Nlcw0KDQpTaSBsYSByZXNwdWVzdGEgcXVlIGJ1c2NhbW9zIGVzIHVuIGZhY3RvciBjb24gbWFzIGRlIDIgdmFsb3JlcyAqc3ZtKCkqIHJlYWxpemFyw6EgdW5hIGNsYXNpZmljYWNpw7NuIG11bHRpY2xhc2UNCg0KYGBge3J9DQpzZXQuc2VlZCAoMSkNCnggPSByYmluZCh4LCBtYXRyaXggKHJub3JtICg1MCAqIDIpICwgbmNvbCA9IDIpKQ0KeSA9IGMoeSwgcmVwICgwICwgNTApKQ0KeFt5ID09IDAgLCAyXSA9IHhbeSA9PSAwICwgMl0gKyAyDQpkYXQgPSBkYXRhLmZyYW1lKHggPSB4LCB5ID0gYXMuZmFjdG9yICh5KSkNCnBhcihtZnJvdyA9IGMoMSwgMSkpDQpwbG90KHgsIGNvbCA9ICh5ICsgMSkpDQoNCnN2bWZpdCA9IHN2bSgNCiAgeX4uLA0KICBkYXRhID0gZGF0ICwNCiAga2VybmVsID0icmFkaWFsIiwNCiAgY29zdCA9IDEwLA0KICBnYW1tYSA9IDENCikNCnBsb3Qoc3ZtZml0ICwgZGF0KQ0KYGBgDQoNCmxhIGZ1bmNpb24gKnN2bSgpKiB0YW1iaWVuIHB1ZWRlIGhhY2VyIG1hcXVpbmFzIGRlIHNvcG9ydGUgZGUgdmVjdG9yZXMgZGUgcmVncmVzaW9uIHNpIGxhIHZhcmlhYmxlIGVzIHVuIG51bWVyaWNhIHkgbm8gdW4gZmFjdG9yLg0KDQojIyMgOS42LjUgQXBwbGljYXRpb24gdG8gR2VuZSBFeHByZXNzaW9uIERhdGENCg0KdXNhcmVtb3MgZWwgc2V0IGRlIGRhdG9zICpLaGFuKiBhY2VyY2EgZGUgNCB0ZWppZG9zIGNhbmNlcmlnZW5vcy4NCg0KYGBge3J9DQpsaWJyYXJ5IChJU0xSKQ0KbmFtZXMoS2hhbikNCg0KZGltKEtoYW4keHRyYWluKQ0KZGltKEtoYW4keHRlc3QpDQpkaW0oS2hhbiR5dHJhaW4pDQpkaW0oS2hhbiR5dGVzdCkNCg0KdGFibGUoS2hhbiR5dHJhaW4gKQ0KdGFibGUoS2hhbiR5dGVzdCApDQpgYGANCg0KdXNhcmVtb3MgdW4gbW9kZWxvIGxpbmVhbCBwb3IgbGEgZ3JhbiBjYW50aWRhZCBkZSBvYnNlcnZhY2lvbmVzIHF1ZSB0ZW5lbW9zDQoNCmBgYHtyfQ0KDQpkYXQ9ZGF0YS5mcmFtZSh4PUtoYW4keHRyYWluICwgeT1hcy5mYWN0b3IoIEtoYW4keXRyYWluICkpDQpvdXQ9c3ZtKHl+LiwgZGF0YT1kYXQgLCBrZXJuZWwgPSJsaW5lYXIiLGNvc3QgPTEwKQ0Kc3VtbWFyeSAob3V0KQ0KDQp0YWJsZShvdXQkZml0dGVkICwgZGF0JHkpDQpgYGANCg0KIyMjbGliZXJhIG1lbW9yaWENCg0KYGBge3J9DQpnYygpDQpgYGANCg0K