{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Game of Life - Genetic Algorithm (Spanish)\n\n* **V1-V2**: Algoritmo Genético - 15 iteraciones.\n* **V3**: Comparación con otros métodos y mejora de predicciones.\n* **V4-V5**: Optimización de los resultados de **V3**."},{"metadata":{},"cell_type":"markdown","source":"## Objetivo\n\nEl objetivo de este notebook es implementar una solución al problema utilizando algoritmos genéticos, de acuerdo al post de [Praveen Kumar](https://www.kaggle.com/c/conways-reverse-game-of-life-2020/discussion/179872).\n\nAcontinuación se muestra un análisis exploratorio intentando relacionar el estado **start** con el estado **stop**, para luego aplicar un algoritmo genético para buscar la mejor solución posible al problema. Luego una comparación con otros métodos y una mejora a la predicción realizada."},{"metadata":{},"cell_type":"markdown","source":"## Análisis exploratorio\n\nLos datos de entrenamiento tienen 3 conjuntos de variables relevantes. La primera es **delta** que representa la iteración desde **start** a **stop**, va de 1 a 5. Luego están los estados **start** y **stop**, cada uno con 625 (25x25) elementos."},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"library(tidyverse)\nlibrary(ggcorrplot)\ninstall.packages(\"genalg\")\nlibrary(genalg)\nlibrary(gridExtra)\nlibrary(ggExtra)\nlibrary(EBImage)\nlibrary(reshape2)\nlibrary(future)\nlibrary(repr)\n\noptions(repr.plot.width=15, repr.plot.height = 10)\n\ntrain <- read_csv(paste0(\"../input/\", list.files(path = \"../input\")[1], \"/train.csv\"))\nhead(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"En el gráfico siguiente se muestra la relación entre la suma de **stop**, la suma de **start** y la diferencia entre ambas por cada **delta**. La diferencia está normalizada por la mediante la suma de ambas variables para evitar la influencia del número de celulas vivas.\n\nSe observa que esta relacion es relativamente lineal y que a medida que se hace más grande **stop**, aumenta la dispersión en torno a la recta que representa a los datos.\n\nAdemás, observando la recta de regresión, se tiene que a mayor **delta** menos cantidad de celulas sobreviven y que el rango de variación es mayor a medida que aumenta **delta**.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"delta <- train %>% select(grep(\"delta\", colnames(train))) %>% as.matrix()\nstart <- train %>% select(grep(\"start_\", colnames(train))) %>% as.matrix()\nstop <- train %>% select(grep(\"stop_\", colnames(train))) %>% as.matrix()\n\ndifSS <- data.frame(delta = factor(delta),\n                    difSS = (start - stop) %>% rowSums(),\n                    stop = stop %>% rowSums(),\n                    start = start %>% rowSums(),\n                    stringsAsFactors = FALSE)\n\ndifSS %>% mutate(difSS = difSS/(stop + start)) %>% \n      ggplot(aes(y = stop, x = start, color = difSS))+\n      geom_point(size = 1) + \n      scale_color_gradient2(low = \"blue\", mid = \"orange\", high = \"blue\",\n                           midpoint = mean(range(difSS$difSS/(difSS$stop + difSS$start))))+\n      geom_smooth(method = \"glm\", color = \"black\", size = 0.5, n = 1000) +\n      facet_wrap(delta~., labeller = label_both) + \n      ggtitle(\"Relación entre Start y Stop\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"En el boxplot se observa un leve aumento en la diferencia cuando aumenta **delta**. Esto hace más complicado estimar una ubicación cuando **delta** crece.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"difSS %>% mutate(difSS = difSS/(stop + start)) %>% \n      mutate(delta = delta %>% as.character()) %>% \n      ggplot(aes(x = delta, y = difSS, fill = delta))+\n      scale_fill_manual(values = colorRampPalette(c(\"blue\",\"orange\",\"green\"))(6))+\n      geom_violin(alpha = 0.6) + \n      geom_point(stat= \"summary\", fun=mean, shape=16, size=1.5, color=\"red\") + \n      ggtitle(\"Boxplot de diferencia entre Start y Stop\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Algoritmo Genético\n\nPara el desarrollo del algoritmo genético se utiliza la librería [genalg](https://cran.r-project.org/web/packages/genalg/index.html). Esta librería tiene la particularidad de que es posible incluir cromosomas iniciales tanto aleatorios, como específicos. Se utiliza esta característica para acotar la zona de busqueda del algoritmo."},{"metadata":{},"cell_type":"markdown","source":"### Zonas de Busqueda\n\nLa tendencia del algoritmo del juego de la vida es a generar modificaciones en los vecinos cercanos, esto es en una vecindad de 3x3 cuando **delta** = 1, hasta llegar a una vecindad máxima de 6x6, por lo que para definir zonas de busque se utilizan las zonas cercanas a las células vivas y así minimizar el tiempo de busqueda de una solución."},{"metadata":{"trusted":true},"cell_type":"code","source":"zoneInv <- function(x, n = 1, umbral = TRUE){\n   mask <- matrix(c(1/8,1/8,1/8,1/8,0,1/8,1/8,1/8,1/8),3,3)\n   # deja guardado el primer filtro\n   x2 <- x %>% filter2(filter =  mask)\n   # itera la convoluvion\n   for (k in c(1:n)) {x <- x %>% filter2(filter =  mask)}\n   x <- sqrt(abs(x*x2)) # suaviza el valor resultante, cuando n muy grande tiende a abarcar todo o nada\n   if(umbral){x <- 1*(x>((0.5)^n))} # umbral para seleccionar la zona\n   return(x)\n}\nzoneInvD <- function(x, n = 2, umbral = T, r = 1){\n   # calcula posiciones de mascaras (combinatoria n entre 8)\n   if(n>0){maskn <- combn(8,n) %>% apply(2,function(x) ifelse(x>=5, x+1, x))}\n   if(n==0){d_maskn <- 1}\n   if(n==1){d_maskn <- 8}\n   if(n>1){d_maskn <- dim(maskn)[2]}\n   # itera la convolucion mascaras n vecinos\n   for (j in 1:d_maskn) {\n      cm <- rep.int(0,9)\n      if(n==1){cm[maskn[j]] <- 1/n}\n      if(n>1){cm[maskn[,j]] <- 1/n}\n      cm <- cm %>% matrix(nrow = 3, ncol = 3)\n      xn <- x %>% filter2(filter = cm)\n      if(j==1){xxn <- xn %>% matrix(nrow = 1)}\n      if(j>1){xxn <- rbind(xxn, xn %>% matrix(nrow = 1))}\n   }\n   xxn <- c(1:d_maskn) %>% sapply(function(z){xxn[z,]*1}) %>% t() %>% colMeans()\n   if(umbral) {xxn <- 1*(xxn>=(r*sum(x)/625))}\n   return(xxn)\n}\n\nset.seed(1)\nt = sample.int(50000,1)\nsugg <- stop[t,] %>% matrix(nrow = 25, ncol = 25) %>% zoneInv(n=1, umbral = F) %>% matrix(nrow = 1)\nsugg <- rbind(sugg, start[t,], stop[t,], stop[t,] - start[t,])\n\n\ngrid.arrange(\n      ggcorrplot(sugg[1,] %>% matrix(nrow = 25, ncol = 25), show.legend = F, title = \"Conv(stop) delta = 1\"),\n      ggcorrplot(sugg[2,] %>% matrix(nrow = 25, ncol = 25), show.legend = F, title = \"Start Original\"),\n      ggcorrplot(sugg[3,] %>% matrix(nrow = 25, ncol = 25), show.legend = F, title = \"Stop Original\"),\n      ggcorrplot(sugg[4,] %>% matrix(nrow = 25, ncol = 25), show.legend = F, title = \"Stop - Start\"),\n      nrow = 2,\n      top = paste0(\"Posibles zonas de busqueda\") \n      )\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Modelo\n\nUna vez definidas las zonas de busqueda se realiza un cálculo de prueba con 50 muestras y 15 iteraciones en promedio (las iteraciones son variables, dependen de **delta** y del número de células vivas de **stop**)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# funcion del game of life hacia adelante\ngameLifeStep <- function(im, n=1){\n      for (k in c(1:n)) {\n            x <- rbind(cbind(im,im,im),\n                       cbind(im,im,im),\n                       cbind(im,im,im))\n            x2 <- matrix(0, ncol = ncol(im), nrow = nrow(im))\n            for (i in (nrow(im)+1):(2*nrow(im))) {\n                  for (j in (ncol(im)+1):(2*ncol(im))) {# sobrepoblacion y subpoblacion por descarte\n                        if(x[i,j]==0){\n                              if(sum(x[(i-1):(i+1), (j-1):(j+1)])==3){\n                                    x2[(i-nrow(im)),(j-ncol(im))] = 1 # reproduccion\n                              }\n                        }\n                        if(x[i,j]==1){\n                              if((sum(x[(i-1):(i+1), (j-1):(j+1)])==3)|(sum(x[(i-1):(i+1), (j-1):(j+1)])==4)){\n                                    x2[(i-nrow(im)),(j-ncol(im))] = 1 # Estasis\n                              }\n                        }\n                  }\n            }\n            im <- x2\n      }\n      return(x2)\n}\n\n# funcion objetivo\nfuncObj <- function(x, y, n=1){\n      x <- x %>% matrix(nrow = 25, ncol = 25) %>% gameLifeStep(n = n) %>% matrix(nrow = 1)\n      return(sum(x!=y)/625)\n}\n# funcion que hace pruebas a un conjunto de stop´s\ntestAG <- function(x, delta, popSizeFactor = 1.25, iters = 10, \n                   mutationChance = 0.005, elitism = NA, verbose = T, suggN = NULL){\n      num = dim(x)[1]\n      vars = dim(x)[2]\n      pop <- matrix(nrow = 1, ncol = vars)\n      evalHist <- matrix(nrow = 1, ncol = vars)\n      eval <- c()\n      model <- list()\n      if(verbose){pb <- txtProgressBar(min = 0, max = num, style = 3)}\n      \n      for (k in 1:num) {\n            z <-  x[k,]\n            sugg <- z\n            sugg <- rbind(sugg, z %>% matrix(nrow = 25, ncol = 25) %>% zoneInvD(n=2, umbral = T)) #2\n            sugg <- rbind(sugg, z %>% matrix(nrow = 25, ncol = 25) %>% zoneInvD(n=3, umbral = T)) #3\n            for (kk in 1:5) {\n               sugg <- rbind(sugg, z %>% matrix(nrow = 25, ncol = 25) %>% zoneInv(n=kk, umbral = T)\n                             %>% matrix(nrow = 1)) # 4-8\n            }\n            if(!is.null(suggN)){sugg <- rbind(sugg, suggN[k,])}\n            \n            pSize <- round(popSizeFactor*dim(sugg)[1])\n            if(pSize<=5){pSize <- 6}\n\n            myProb <-rbga.bin(size=vars, \n                              evalFunc = function(x, n=delta[k]){funcObj(x, y = z, n = n)}, \n                              showSettings = F, \n                              popSize = pSize, \n                              iters = ceiling(1.5*sqrt(iters)*sqrt(delta[k]+iters)*log1p(sum(z)*iters)),\n                              suggestions = sugg,\n                              mutationChance = mutationChance,\n                              elitism = elitism,\n                              verbose = F)\n            \n            model[[k]] <- myProb\n            pop <-rbind(pop, \n                        myProb$population[order(myProb$evaluations)[1],])\n            \n            if(k<=1){evalHist <- matrix(myProb$evaluations, nrow = 1)}\n            if(k>1){evalHist <- rbind(evalHist, myProb$evaluations)}\n            \n            eval <- c(eval, \n                      funcObj(x = matrix(myProb$population[order(myProb$evaluations)[1],], nrow = 25, ncol = 25),\n                              y = z,\n                              n = delta[k]))\n            if(verbose){setTxtProgressBar(pb, k)}\n      }\n      pop <- pop[-1,]\n      if(verbose){close(pb)}\n      return(list(eval = eval, pop = pop, evalHist = evalHist, model = model))\n}\n\nset.seed(2)\nid <- sample.int(50000,50)\np <- testAG(x = stop[id,], delta = delta[id], popSizeFactor = 1.1,\n            mutationChance = 0.03, iters = 5*3, elitism = T, verbose = F)\nsummary(p$eval)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data.frame(loss = p$eval, delta = factor(delta[id])) %>% \n   ggplot(aes(y = loss, x = delta, fill = delta)) +\n   geom_violin(alpha = 0.6)+ ylab(\"loss (1-accuracy)\")+\n   scale_fill_manual(values = colorRampPalette(c(\"blue\",\"orange\",\"green\"))(6))+\n   geom_boxplot(width = 0.1, outlier.shape = NA) +\n   ggtitle(\"Boxplot del error en función de delta\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"i1 = which(p$eval==min(p$eval))[1]\ni2 = which(p$eval==max(p$eval))[1]\ni3 = which(p$eval==p$eval[order(p$eval)][round(0.5*length(p$eval))])[1]\ngrid.arrange(\n      ggcorrplot(p$pop[i1,] %>% matrix(nrow = 25, ncol = 25) %>% \n                       gameLifeStep(n = delta[id[i1]]), \n                 show.legend = F, \n                 title = paste0(\"Prediccion (min)\")),\n      ggcorrplot(unlist(stop[id[i1],]) %>% matrix(nrow = 25, ncol = 25), \n                 show.legend = F, title = paste0(\"Original\")),\n      data.frame(mean = p$model[[i1]]$mean, \n                 best = p$model[[i1]]$best, \n                 generation = c(1:length(p$model[[i1]]$mean))) %>% \n         melt(id = \"generation\") %>% \n         ggplot(aes(y = value, x = generation, colour = variable)) + \n         geom_point() + ylab(\"loss\") + ggtitle(paste0(\"Loss: \", round(min(p$eval), 2))) + theme_minimal(),\n      \n      \n      ggcorrplot(p$pop[i2,] %>% matrix(nrow = 25, ncol = 25) %>% \n                       gameLifeStep(n = delta[id[i2]]), \n                 show.legend = F, \n                 title = paste0(\"Prediccion (max)\")),\n      ggcorrplot(unlist(stop[id[i2],]) %>% matrix(nrow = 25, ncol = 25), \n                 show.legend = F, title = paste0(\"Original\")),\n      data.frame(mean = p$model[[i2]]$mean, \n                 best = p$model[[i2]]$best, \n                 generation = c(1:length(p$model[[i2]]$mean))) %>% \n         melt(id = \"generation\") %>% \n         ggplot(aes(y = value, x = generation, colour = variable)) + \n         geom_point() + ylab(\"loss\") + ggtitle(paste0(\"Loss: \", round(max(p$eval), 2))) + theme_minimal(),\n      \n      \n      ggcorrplot(p$pop[i3,] %>% matrix(nrow = 25, ncol = 25) %>% \n                       gameLifeStep(n = delta[id[i3]]), \n                 show.legend = F, \n                 title = paste0(\"Prediccion (med)\")),\n      ggcorrplot(unlist(stop[id[i3],]) %>% matrix(nrow = 25, ncol = 25), \n                 show.legend = F, title = paste0(\"Original\")),\n      data.frame(mean = p$model[[i3]]$mean, \n                 best = p$model[[i3]]$best, \n                 generation = c(1:length(p$model[[i3]]$mean))) %>% \n         melt(id = \"generation\") %>% \n         ggplot(aes(y = value, x = generation, colour = variable)) + \n         geom_point() + ylab(\"loss\") + \n         ggtitle(paste0(\"Loss: \", round(p$eval[order(p$eval)][round(0.5*length(p$eval))], 2))) + theme_minimal(),\n      \n      \n      nrow = 3,\n      top = \"Visualizacion de Stop de valores extremos (minimmo, maximo, mediano)\" \n      )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Modelo aplicado a test\n\nPara los cálculos de **test** se utiliza la libreria [future](https://cran.r-project.org/web/packages/future/vignettes/future-1-overview.html) para paralelizar los cálculos. Considerando 15 iteraciones la optimización completa demora alrededor de 15 hrs."},{"metadata":{"trusted":true},"cell_type":"code","source":"test <- read_csv(paste0(\"../input/\", list.files(path = \"../input\")[1], \"/test.csv\"))\n\ndelta <- test %>% select(grep(\"delta\", colnames(test)))\nstop <- test %>% select(grep(\"stop_\", colnames(test)))\n \n# if(!dir.exists(\"AG test 2/\")){dir.create(\"AG test 2/\")}\n# \n# it <- 5*3\n# \n# seqk <- seq.int(131,250,5)\n# for (k in seqk) {# debe demorar 15 hrs aprox\n#    if(k==seqk[1]){cat(\"================= START =================\\n\")}\n#    t1 <- Sys.time()\n# \n#    id1 <- c(1:200) + 200*(k-1) # debe demorar alrededor de 5 mins\n#    id2 <- c(1:200) + 200*(k) # debe demorar alrededor de 5 mins\n#    id3 <- c(1:200) + 200*(k+1) # debe demorar alrededor de 5 mins\n#    id4 <- c(1:200) + 200*(k+2) # debe demorar alrededor de 5 mins\n#    id5 <- c(1:200) + 200*(k+3) # debe demorar alrededor de 5 mins\n#    \n#    p1 <- future({ testAG(x = stop[id1,], delta = delta$delta[id1], popSizeFactor = 1.25,\n#                          mutationChance = 0.03, iters = it, verbose = F) }) %plan% multiprocess\n#    p2 <- future({ testAG(x = stop[id2,], delta = delta$delta[id2], popSizeFactor = 1.25,\n#                          mutationChance = 0.03, iters = it, verbose = F) }) %plan% multiprocess\n#    p3 <- future({ testAG(x = stop[id3,], delta = delta$delta[id3], popSizeFactor = 1.25,\n#                          mutationChance = 0.03, iters = it, verbose = F) }) %plan% multiprocess\n#    p4 <- future({ testAG(x = stop[id4,], delta = delta$delta[id4], popSizeFactor = 1.25,\n#                          mutationChance = 0.03, iters = it, verbose = F) }) %plan% multiprocess\n#    p5 <- future({ testAG(x = stop[id5,], delta = delta$delta[id5], popSizeFactor = 1.25,\n#                          mutationChance = 0.03, iters = it, verbose = F) }) %plan% multiprocess\n# \n#    pp1 <- value(p1)\n#    pp2 <- value(p2)\n#    pp3 <- value(p3)\n#    pp4 <- value(p4)\n#    pp5 <- value(p5)\n# \n#    saveRDS(pp1, paste0(\"AG test 2/AG_test_\", k, \".rds\"))\n#    saveRDS(pp2, paste0(\"AG test 2/AG_test_\", k+1, \".rds\"))\n#    saveRDS(pp3, paste0(\"AG test 2/AG_test_\", k+2, \".rds\"))\n#    saveRDS(pp4, paste0(\"AG test 2/AG_test_\", k+3, \".rds\"))\n#    saveRDS(pp5, paste0(\"AG test 2/AG_test_\", k+4, \".rds\"))\n# \n#    t2 <- Sys.time()\n#    \n#    cat(paste0(\"---- iter \", k,\" ----\", \"\\n\"))\n#    print(t2 -t1)\n#    cat(\"                                \\n\")\n#    cat(paste0(\"loss \",k,\": \", summary(pp1$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",k+1,\": \", summary(pp2$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",k+2,\": \", summary(pp3$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",k+3,\": \", summary(pp4$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",k+4,\": \", summary(pp5$eval)[4]), \"\\n\")\n#    cat(\"=========================================\\n\")\n# \n# }\npath <- \"../input/conways-reverse-game-of-life-2020-ag1/AG test 2/AG_test_\"\n\neval <- c()\nfor (i in 1:250) {eval <- c(eval, readRDS(paste0(path, i, \".rds\"))$eval)}\nsummary(eval, digits = 6)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"El error promedio anda en torno al 10%. Se observa que el error aumenta a medida que crece **delta**, por lo que es necesario darle más tiempo de cálculo a los delta más grandes."},{"metadata":{"trusted":true},"cell_type":"code","source":"eval <- data.frame(delta = factor(delta$delta), eval = eval)\n\neval %>% \n   ggplot(aes(y = eval, x = delta, fill = delta)) +\n   geom_violin(alpha = 0.4)+ ylab(\"loss (1-accuracy)\")+\n   scale_fill_manual(values = colorRampPalette(c(\"blue\",\"orange\",\"green\"))(6))+\n   geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.7) +\n   ggtitle(\"Distribución del error en función de delta (test)\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n### Comparación con otros modelos.\n\nPara la comparación considero modelos distintos a los algoritmos genéticos, particularmente se consideran los modelos de los siguientes notebooks:\n\n* Game of Life - Genetic Algorithm (Spanish) (**V-2**) - error: 0.10234.\n* [Game of Life - Z3 Constraint Satisfaction](https://www.kaggle.com/c/conways-reverse-game-of-life-2020/notebooks?sortBy=scoreAscending&group=everyone&pageSize=20&competitionId=22122) (**V-31**) - error: 0.08679.\n* [[CRGL] Probability Extension | True Target Problem](https://www.kaggle.com/yakuben/crgl-probability-extension-true-target-problem) (**V-19**) - error: 0.11117.\n\nSeleccioné los 2 modelos anteriores por ser muy distintos a los algoritmos genéticos y porque me pareció interesante el enfoque de cada uno.\n\nCon estos 3 modelos se seleccionarán las predicciones de **start** que tengan el menor error, el resultado se muestra a continuación."},{"metadata":{"trusted":true},"cell_type":"code","source":"loss <- function(x, y = unlist(stop[6,]), n=1){\n      x <- x %>% matrix(nrow = 25, ncol = 25) %>% gameLifeStep(n = n) %>% matrix(nrow = 1)\n      return(sum(x!=y)/625)\n}\n\ntest <- read_csv(paste0(\"../input/\", list.files(path = \"../input\")[1], \"/test.csv\"))\n\ndelta <- test %>% select(grep(\"delta\", colnames(test))) %>% as.matrix()\nstop <- test %>% select(grep(\"stop_\", colnames(test))) %>% as.matrix()\n\npath <- \"../input/conways-reverse-game-of-life-2020-ag1/\"\n\nsub_desareca <- (read_csv(paste0(path, \"sub_desareca.csv\")) %>% as.matrix())[,-1]\nsub_McGuigan_Z3 <- (read_csv(paste0(path, \"sub_McGuigan_Z3.csv\")) %>% as.matrix())[,-1]\nsub_yakuben_prob <- (read_csv(paste0(path, \"sub_yakuben_prob.csv\")) %>% as.matrix())[,-1]\n\neval_desareca <- c(1:50000) %>% sapply(function(z){loss(x = sub_desareca[z,], y = stop[z,], n = delta[z])})\neval_McGuigan_Z3 <- c(1:50000) %>% sapply(function(z){loss(x = sub_McGuigan_Z3[z,], y = stop[z,], n = delta[z])})\neval_yakuben_prob <- c(1:50000) %>% sapply(function(z){loss(x = sub_yakuben_prob[z,], y = stop[z,], n = delta[z])})\n\nsubComp <- data.frame(samp = c(1:50000), delta = delta,\n                      Des = eval_desareca,\n                      Mcg = eval_McGuigan_Z3,\n                      Yak = eval_yakuben_prob,\n                      stringsAsFactors = F)\n\n# Identifica cual modelo genera una mejor aproximacion\nidEval <- subComp %>% select(Des, Mcg, Yak) %>% apply(1, which.min)\n\n#trae el menor error entre los modelos.\neval <- c(1:50000) %>% sapply(function(x){(subComp %>% select(Des, Mcg, Yak))[x, idEval[x]]})\n\neval %>% summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"El resultado es considerablemente mejor que cualquiera de los 3 modelos por separado.\n\nA continuación se muestra una comparación, en ella se observa que en general los 3 métodos abordan condiciones distintas, tanto para promedio de células vivas (**meanStop**) como para **delta**. Esta comparación muestra que los 3 métodos se pueden complementar muy bien, lo que concuerda con la mejora mostrada al fusionar los resultados. Los mejores resultados están mayoritariamente entre el algortimo genético y el Z3 solver."},{"metadata":{"trusted":true},"cell_type":"code","source":"subComp <- subComp %>% mutate(minEval = c(\"Des\", \"Mcg\", \"Yak\")[idEval])\nsubComp <- subComp %>% mutate(meanStop = stop %>% rowMeans())\n\n(subComp %>% \n      mutate(minEval = factor(minEval, labels = c(\"Genetic Algorithm\", \"Z3 Solver\", \"Prob. Extension\"))) %>% \n      select(minEval) %>% table())/dim(subComp)[1]\n\nsubComp %>% mutate(minEval = factor(minEval, labels = c(\"Genetic Algorithm\", \"Z3 Solver\", \"Probability Extension\"))) %>% \n      ggplot(aes(y = meanStop, minEval, fill = minEval)) + \n      geom_violin(alpha = 0.4) + \n      geom_boxplot(width = 0.1, outlier.shape = T, alpha = 0.6) +\n      scale_fill_manual(values = colorRampPalette(c(\"blue\",\"orange\",\"green\"))(3))+\n      facet_wrap(delta~., labeller = label_both)+\n      ggtitle(\"Comparación entre distintos métodos\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Optimización de resultados comparados\n\nAhora nuevamente optimizamos los últimos resultados con un algoritmo genético. En este caso se considera optimizar sólo las 22000 muestras que tienen un error mayor a 0 para disminuir los tiempos de cálculo.\n\nSe consideran 50 iteraciones promedio por cada muestra, un factor de mutación de 0.01 y cromosoma inicial a los resultados anteriores. El tiempo de ejecución de este algoritmo es de 22 hrs aproximadamente. \n\nSe observa una mejora en los resultados, pasando de 0.06590 a 0.06189."},{"metadata":{"trusted":true},"cell_type":"code","source":"path <- \"../input/conways-reverse-game-of-life-2020-ag1/\"\n# ordena los indices por error\n# subId <- data.frame(id=c(1:50000), eval = eval) %>% arrange(-eval)\nsubId <- readRDS(paste0(path, \"subId_EvalMin.rds\"))\n\n# selecciona la prediccion con menor error\nstartPredMin <- c(1:50000) %>% \n   sapply(function(x){\n      if(idEval[x]==1){y <- sub_desareca[x,]}\n      if(idEval[x]==2){y <- sub_McGuigan_Z3[x,]}\n      if(idEval[x]==3){y <- sub_yakuben_prob[x,]}\n      return(y)\n   }) %>% t()\n\n# Algoritmo genetico para 22000 muestras con error mayor a 0\n\n# if(!dir.exists(\"AG test 3/\")){dir.create(\"AG test 3/\")}\n# it <- 50\n# seqk <- subId$id\n# for (k in 16:22) {# quedé por comenzar iter:8 # hasta 22 iteraciones, debe demorar 22.5 hrs aprox\n#    if(k==8){cat(\"================= START =================\\n\")}\n#    t1 <- Sys.time()\n#    \n#    id1 <- seqk[c(1:200) + 1000*(k-1)] \n#    id2 <- seqk[c(201:400) + 1000*(k-1)]\n#    id3 <- seqk[c(401:600) + 1000*(k-1)]\n#    id4 <- seqk[c(601:800) + 1000*(k-1)]\n#    id5 <- seqk[c(801:1000) + 1000*(k-1)]\n#    \n#    p1 <- future({ testAG(x = stop[id1,], delta = delta[id1], popSizeFactor = 1.1,\n#                          mutationChance = 0.01, iters = it, verbose = F,\n#                          suggN = startPredMin[id1,]) }) %plan% multiprocess\n#    p2 <- future({ testAG(x = stop[id2,], delta = delta[id2], popSizeFactor = 1.1,\n#                          mutationChance = 0.01, iters = it, verbose = F,\n#                          suggN = startPredMin[id2,]) }) %plan% multiprocess\n#    p3 <- future({ testAG(x = stop[id3,], delta = delta[id3], popSizeFactor = 1.1,\n#                          mutationChance = 0.01, iters = it, verbose = F,\n#                          suggN = startPredMin[id3,]) }) %plan% multiprocess\n#    p4 <- future({ testAG(x = stop[id4,], delta = delta[id4], popSizeFactor = 1.1,\n#                          mutationChance = 0.01, iters = it, verbose = F,\n#                          suggN = startPredMin[id4,]) }) %plan% multiprocess\n#    p5 <- future({ testAG(x = stop[id5,], delta = delta[id5], popSizeFactor = 1.1,\n#                          mutationChance = 0.01, iters = it, verbose = F,\n#                          suggN = startPredMin[id5,]) }) %plan% multiprocess\n#    \n#    pp1 <- value(p1)\n#    pp2 <- value(p2)\n#    pp3 <- value(p3)\n#    pp4 <- value(p4)\n#    pp5 <- value(p5)\n#    \n#    saveRDS(pp1, paste0(\"AG test 3/AG_test_\", 5*(k-1)+1, \".rds\"))\n#    saveRDS(pp2, paste0(\"AG test 3/AG_test_\", 5*(k-1)+2, \".rds\"))\n#    saveRDS(pp3, paste0(\"AG test 3/AG_test_\", 5*(k-1)+3, \".rds\"))\n#    saveRDS(pp4, paste0(\"AG test 3/AG_test_\", 5*(k-1)+4, \".rds\"))\n#    saveRDS(pp5, paste0(\"AG test 3/AG_test_\", 5*(k-1)+5, \".rds\"))\n#    \n#    t2 <- Sys.time()\n#    \n#    cat(paste0(\"---- iter \", k,\" ----\", \"\\n\"))\n#    print(t2 -t1)\n#    cat(\"                                \\n\")\n#    cat(paste0(\"loss \",5*(k-1)+1,\": \", summary(pp1$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",5*(k-1)+2,\": \", summary(pp2$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",5*(k-1)+3,\": \", summary(pp3$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",5*(k-1)+4,\": \", summary(pp4$eval)[4]), \"\\n\")\n#    cat(paste0(\"loss \",5*(k-1)+5,\": \", summary(pp5$eval)[4]), \"\\n\")\n#    cat(\"=========================================\\n\")  \n# }\n\n\n# carga los resultados de la optimizacion\n\npath <- \"../input/conways-reverse-game-of-life-2020-ag1/AG test 3/AG_test_\"\n\neval_Next <- c()\nfor (i in 1:110) {eval_Next <- c(eval_Next, readRDS(paste0(path, i, \".rds\"))$eval)}\n\nAG <- matrix(nrow = 0, ncol = 625)\nfor (i in 1:110) {AG <- rbind(AG, readRDS(paste0(path, i, \".rds\"))$pop)}\n\n# Modifica solo los mejores resultados entre la version anterior y la optimización actual.\nsub_comparacion <- read_csv(\"../input/conways-reverse-game-of-life-2020-ag1/sub_comp.csv\") %>% as.matrix()\nstartPred_Next <- sub_comparacion\nfor (i in 1:22000) {\n   k <- subId$id[i]\n   if(eval_Next[i]<subId$eval[i]){startPred_Next[k,-1] <- AG[i,]}\n}\n\neval_x <- c(1:50000) %>% sapply(function(z){loss(x = sub_comparacion[z,-1], y = stop[z,], n = delta[z])})\neval_y <- c(1:50000) %>% sapply(function(z){loss(x = startPred_Next[z,-1], y = stop[z,], n = delta[z])})\n\nrbind(\n   StartBefore = eval_x %>% summary(),\n   StartAfter = eval_y %>% summary()\n   )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"write_csv(data.frame(startPred_Next), \"submission.csv\")","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"3.6.3"}},"nbformat":4,"nbformat_minor":4}