In papabloblo/sudoku: SUDOKU - Práctica 1 de Computación Evolutiva
library(sudoku)
library(tidyverse)
Introducción
Los algoritmos propios de la computación evolutiva pueden aplicarse en campos muy diversos. El objetivo de esta actividad es el de describir el comportamiento de los algoritmos genéticos aplicados a la resolución del pasatiempo conocido como Sudoku.
Comenzaremos describiendo en qué consiste el problema del sudoku y definiéndolo de una manera formal. Después, se comentarácómo se ha llevado a cabo la implementación en R para pasar a hacer una evaluación experimental que nos permita discutir el comportamiento.
Descripción del problema del sudoku
El juego del Sudoku comenzó a ser popular en Japón en 1986 y consiguió popularidad internacional en 2005 [@Lynce2006]. Se trata de un pasatiempo cuyas reglas son simples de entender pero su resolución puede ser todo un reto. Además, un sudoku bien construido tiene una única solución y puede ser resuelto mediante el razonamiento sin tener que recurrir a procedimientos de ensayo y error.
Los sudokus suelen tener asociado un nivel de dificultad que suele depender del número de casillas rellenas en el inicio. En [@McGuire2014] se muestra que el número mínimo de casillas que deben de estar rellenas para que un sudoku tenga solución (única) es 17.
¿En qué consiste?
El juego del Sudoku tradicional se representa por una matriz de dimensión $9x9$ en la que se disponen números enteros del 1 al 9. En un inicio solo algunas de las casillas están ocupadas con números siendo el objetivo del juego rellenar las casillas vacías.
Los elementos principales del sudoku serán las filas, las columnas y las subcuadrículas como se muestra en (imagen).
Para que el sudoku se considere resuelto, se ha de verificar que:
- cada fila contiene una permutación de los 9 enteros apareciendo cada uno de ellos una y solo una vez,
- cada columna contiene una permutación de los 9 enteros apareciendo cada uno de ellos una y solo una vez,
- cada subcuadrícula debe contener una permutación de los 9 enteros apareciendo cada uno de ellos una y solo una vez.
Estas reglas de resolución explican porqué Sudoku significa "un solo número" en japonés [@Lynce2006].
Formalmente
El problema puede ser expresado matemáticamente de múltiples formas. El desarrollo elegido tiene la motivación de que la aplicación del algoritmo evolutivo sea más natural.
El sudoku se se expresa como un vector $\mathbf{s} = (s_1, s_2, \ldots, s_{81})$ donde $s_i \in {0, 1, 2,\ldots,9}$ para $i = 1, \ldots, 81$. Los elementos sobre los que se deben evaluar la correcta resolución del problema y que servirán para calcular la función fitness del algoritmo evolutivo se describen matemáticamente como
$$
columna(j) = {i \in {1,\ldots, 81}\, | \, mod(i, 9)= mod(j, 9)}
$$
$$
fila(j) = {i \in {1,\ldots, 81}\, | \, \lfloor j/9 \rfloor = \lfloor i/9\rfloor}
$$
donde $\lfloor a \rfloor$ se refiere a la parte entera de $a$.
Por otro lado, definimos la funciones:
$$
c(i) = \sum_{{j \in columna(i)|j \ne i}} I(s_i = s_j)
$$
$$
f(i) = \sum_{{j \in fila(i)|j \ne i}} I(s_i = s_j) \
$$
$$
s(i) = \sum_{{j \in subcuadrícula(i)|j \ne i}} I(s_i = s_j) \
$$
El sudoku se considerará resuelto si y solo si $\sum_i f(i) = 0$, $\sum_i c(i) = 0$ y $\sum_i s(i) = 0$, o lo que es lo mismo
$$
\sum_{i=1}^{81} \left[f(i) + c(i) + s(i)\right] = 0
$$
Podemos definir la siguiente función
$$
fitness(s) = \sum_{i=1}^{81} f(i) + c(i) + s(i)
$$
Cualquier configuración de números que no resuelva el sudoku tendrá un $$fitness(s) > 0$$ por lo que, la solución que resuelve el sudoku se puede expresar como
$$
s^* = \underset{s}{\arg\min} \, fitness(s)
$$
-->
Implementación
El juego del sudoku se va a abordar mediante un algoritmo genético. Hay distintas estrategias que se podrían seguir para este problema. En este caso, la implementación se ha realizado conforme a las características que se detallan a continuación.
Estructura de códigos en R
Los códigos implementados se encuentran en la carpeta \R.
Para que los resultados de las siguientes secciones sea reproducibles, se fija la semilla para la generación de valores aleatorios
set.seed(123)
Representación de cada individuo
Cada celda $s_i$ del sudoku estará asociada a un gen representado por números enteros del 1 al 9, es decir, $s_i\in {1, 2, \ldots, 9}$. El genotipo estará entonces representado por una cadena de enteros. Lo que es lo mismo, se realizará una muestra con reemplazamiento de los enteros del 1 al 9 con una longitud determinada.
La función random_integer_representation() es la que permite obtener este tipo de representación. Necesita los argumentos:
valores: vector con los valores posibles que puede tomar cada gen.tam: número entero con el tamaño del genotipo (número de genes).
Como ejemplo:
random_integer_representation(valores = 1:9, tam = 9)
Aunque el genotipo de cada individuo contendrá un total de $81$ genes, solo se aplicará el algoritmo sobre aquellos genes (casillas) que estén vacíos en la solución inicial. Por tanto, el algoritmo genético manejará el genotipo $s_0$ de forma que
$$
s_i = s^1_i si
$$
Función de adaptación (fitness)
La función de fitness que el algoritmo tratará de minimizar será la que aparece en la ecuación ]bla]. Como se dijo anteriormente, un valor de $fitness(s)=0$ asegura que el juego ha sido resuelto. El valor máximo de esta función se alcanzará cuando todas las casillas (genes) contengan el mismo valor y tendrá un valor de $fitness(s) = 972$
Las funciones implementadas para el cálculo de la función de fitness se encuentran en \R\fitness.R.
La función fitness_sudoku() calcula la función de fitness para un individuo dado y necesita los argumentos:
x: un vector con el genotipo del individuo.ind_cuadricula: un vector indicando qué posiciones equivalen a la representación de las subcuadrículas.solucion_inicial: un vector representando la solución inicial del sudoku que contendrá NA en aquellas posiciones que originalmente aparecían vacías.
En nuestro caso, el ind_cuadricula se representa como
ind_cuadricula <- c(
rep(rep(x = 1:3, each = 3), 3),
rep(rep(x = 4:6, each = 3), 3),
rep(rep(x = 7:9, each = 3), 3))
ind_cuadricula
Viendo este vector como una matriz, se pone de relieve cómo están representadas las subcuadrículas
matrix(ind_cuadricula, nrow = 9, ncol = 9)
Un ejemplo de solución inicial del sudoku se podría representar así:
sudoku_facil <- c(NA, NA, NA, NA, 3, NA, NA, NA, 4,
NA, 9, NA, 4, NA, 6, NA, 7, NA,
NA, 5, NA, NA, NA, NA, 3, 8, NA,
NA, NA, NA, NA, 7, 8, NA, NA, 3,
3, NA, NA, NA, NA, NA, 6, 9, NA,
5, 4, NA, 6, NA, NA, NA, 2, NA,
7, NA, 5, NA, 2, 4, NA, NA, NA,
9, 8, 4, NA, 6, 5, 2, NA, NA,
NA, 2, 6, NA, 8, NA, NA, NA, 9)
El genotipo de los individuos estará representado por el número de casillas vacías en la solución inicial, es decir,
tam_individuos <- sum(is.na(sudoku_facil))
tam_individuos
Simulamos un individuo con 48 genes, en este caso,
individuo <- random_integer_representation(valores = 1:9, tam = tam_individuos)
Ya podemos calcular el valor de la función de fitness para este individuo.
fitness_sudoku(x = individuo,
ind_cuadricula = ind_cuadricula,
solucion_inicial = sudoku_facil)
Inicialización de la población
La población inicial del algoritmo se realizatá de forma aleatoria, es decir, cada individuo estará representado por los genes que le corresponda y cada gen se inicializará como un número entero aleatorio entre 1 y 9.
Esto se permite gracias a la función generacion_poblacion() cuyos argumentos son:
valores_posibles: vector con los valores posibles que puede tomar cada gen.num_genes: número entero con el tamaño del genotipo (número de genes).tam_poblacion: entero con el número de individuos de los que estará formada la población.
El resultado es una lista compuesta por los distintos individuos. Por ejemplo,
poblacion <- generacion_poblacion(valores_posibles = 1:9,
num_genes = tam_individuos,
tam_poblacion = 3)
poblacion
Selección de padres
Para poder realizar el cruce entre individuos que se describe en la siguiente sección, es necesario elegir los padres, es decir, aquellos individuos entre los que realizar el cruce. En este caso se realiza la selección por torneo.
Tenemos la función seleccion_padres() cuyos argumentos son:
num_padres:poblacion: lista con los individuos que compoennfitness_poblacion: vector con el valor de fitness para cada genotipo de la población.k: entero con el número de competidores del torneo
El resultado es una lista con los padres elegidos por torneo.
padres <- seleccion_padres(num_padres = 2,
poblacion = poblacion,
fitness_poblacion = sapply(poblacion,
fitness_sudoku,
ind_cuadricula = ind_cuadricula,
solucion_inicial = sudoku_facil),
k = 2)
padres
Cruce
El cruce entre los padres se realizará mediante cruce por un punto.
Función one_point_crossover() con argumentos:
padres: lista con los padres seleccionados para el cruce.prob_cruce: numeric en [0, 1] representando la probabilidad de que los padres se crucen y tengan descendencia o se mantengan.
resultado_cruce <- one_point_crossover(padres = padres, prob_cruce = 0.75)
resultado_cruce
Mutación
Una vez producido el cruce entre padres y obtenida su descendencia, se aplicará mutación random resetting, es decir, cada uno de los genes de cada individuo podrá sufrir una mutación con una probabilidad dada. En el caso de que se produzca una mutación en un gen, esta podrá tomar un valor entero aleatorio entre 1 y 9.
random_resetting() con parámetros:
x: individuo sobre el que provocar mutacionesprob: probabilidad de que se aplique mutación en cada gen.valores_posibles: vector con valores posibles para la mutación de cada gen.
random_resetting(individuo, prob = 0.1, valores_posibles = 1:9)
Selección de supervivientes
Se utilizará un modelo de tipo generacional, reemplazando la población actual enteramente por una nueva aplicando, eso sí, elitismo para conservar el emjor individuo en la siguiente generación y eliminando el peor.
El algoritmo
Una ejecución del algoritmo se realizará gracias a la función algoritmo_genetico() con parámetros
poblacion_inicial: lista con los individuos que compone la población inicial.funcion_fitness: función que calculael fitness de cada individuo de la población.genes_fijos: ...valores_mutacion: vector con los valores posibles para la mutación.num_padres: entero con el número de padres.prob_cruce: probabilidad de que se produzca cruce entre padres.tam_torneo: tamaño del torneo en la selección de padres por torneo.prob_mutacion: probabilidad de que se produzca una mutación en cada gen del individuo.tam_poblacion: tamaño de la población ¿?max_iter: iteraciones máximas que realizará el algoritmoprint_each: imprimir por pantalla información del algoritmo cada print_each iteraciones.
Se ha elegido como lenguaje de programación R. R es un tipo de lenguaje Open Source con una gran aceptación dentro del área del data science.
Todas las funciones se han implementado de tal forma que puedan utilizarse para aplicar un algoritmo genético en general y no se han particularizado para el problema del Sudoku a excepción de la función de fitness.
Toda la documentación de cada una de las funciones implementadas se encuentra en el pdf adjunto.
Evaluación experimental
Los algoritmos genéticos llevan asociados una serie de parámetros que rigen su comportamiento y pueden impactar en los resultados obtenidos. Por lo tanto, es fundamental conocer cómo afectan los distintos parámetros al algoritmo para poder tener un conocimiento sobre los mismos que nos permita ser conscientes de los cambios. Por ejemplo, en función del problema, nos puede interesar favorecer una combinación de parámetros que necesite un menor tiempo computacional (un menor número de iteraciones) para obtener una solución cuya función de fitness sea ligeramente peor (en media) que la de otra combinación de parámetros pero que necesite un alto número de iteraciones.
Metodología
Para contrastar el comportamiento del algoritmo genético aplicado al problema de resolver un Sudoku, se han elegido dos configuraciones iniciales, una fácil con 33 casillas iniciales rellenas y otra difícil, con 30 casillas rellenas. Como se dijo en la sección BLA, se considera que la dificultad del problema del Sudoku es inversamente proporcional al número de celdas rellenas al principio.
En la implementación descrita en la sección BLA del algoritmo genético, los parámetros del algoritmo serían:
- Tamaño de la población.
- En la selección de padres por torneo, el tamaño del torneo.
- Probabilidad de cruce
- Probabilidad de mutación.
En nuestro caso, una de las preguntas fundamentales que queremos responder es cómo impactan los diferentes valores de la probabilidad de mutación en el algoritmo genético.
A continuación se describen las dos metodologías seguidas.
Aproximación 1
Se han fijado todos los parámetros del algoritmo a excepción de la probabilidad de mutación para la cual se han elegido valores. Se ha ejecutado el algoritmo 50 veces.
Aproximación 2
Para tratar de obtener un resultado del impacto global de los parámetros, se han generado $50$ configuraciones aleatorias de los parámetros con la siguiente lógica:
- El tamaño de la población se fija en $10$ para la instancia sencilla y $100$ para la instancia compleja.
- Probabilidad de cruce se supone una variable aleatoria uniforme $U(0.5, 1)$.
- Tamaño del torneo: valor entero aleatorio entre $2$ y $10$, ambos incluidos.
- Probabilidad de cruce: valor con distribución uniforme $U(0.5, 1)$.
- Probabilidad de mutación: permutación aleatoria de 50 números igualmente espaciados entre 0 y 0.25.
- $10$ ejecuciones para cada combinación de parámetros
Métricas y gráficas
Instancia sencilla
La instancia sencilla se ha tomado de http://www.sudoku-online.org/ en la categoría fácil 10 (#1217). Su vector se implementa de la siguiente forma
sudoku_facil <- c(NA, NA, NA, NA, 3, NA, NA, NA, 4,
NA, 9, NA, 4, NA, 6, NA, 7, NA,
NA, 5, NA, NA, NA, NA, 3, 8, NA,
NA, NA, NA, NA, 7, 8, NA, NA, 3,
3, NA, NA, NA, NA, NA, 6, 9, NA,
5, 4, NA, 6, NA, NA, NA, 2, NA,
7, NA, 5, NA, 2, 4, NA, NA, NA,
9, 8, 4, NA, 6, 5, 2, NA, NA,
NA, 2, 6, NA, 8, NA, NA, NA, 9)
Nota: tanto la instancia compleja, como la instancia sencilla se recogen el script \sudoku_prueba.R en los objetos sudoku_facil y sudoku_dificil, respectivamente.
Aproximación 1
Se genera un grid con la siguiente configuración de parámetros
grid <- data.frame(tam_poblacion = 10,
prob_cruce = 0.9,
tam_torneo = 2,
prob_mutacion = c(0, 0.01, 0.05, 0.1, 0.15, 0.25))
grid
Como ya se comentó, para esta primera aproximación se mantienen todos los parámetros constantes a excepción de la probabilidad de mutación que toma valores $0$, $0.01$, $0.05$, $0.1$, $0.15$ y $0.25$.
Para cada combinación de parámetros se ha repetido la ejecución $50$ veces dejando un máximo de $5000$ iteraciones.
Esta ejecución se recoge en
library(parallel)
set.seed(1234)
ini <- Sys.time()
resultados_dificil <- mcmapply(pruebas_ga,
num_pruebas = 50,
num_padres = 2,
tam_poblacion = grid$tam_poblacion,
prob_cruce = grid$prob_cruce,
tam_torneo = grid$tam_torneo,
prob_mutacion = grid$prob_mutacion,
max_iter = 5000,
print_each = 250,
MoreArgs = list(generacion_poblacion_ini = generacion_poblacion,
genes_fijos = sudoku_facil,
funcion_fitness = fitness_sudoku,
valores_posibles = 1:9,
valores_mutacion = 1:9
)
,
SIMPLIFY = FALSE
)
fin <- Sys.time()
fin-ini
save(resultados_dificil, file = paste0("data/prueba_facil_", Sys.Date(), ".RData"))
Los resultados de esta ejecución se encuentran en el archivo actividad/prueba_facil_2017-12-10.RData. El objeto resultados_facil contiene una gran cantidad de información (población de cada iteración, fitness de cada individuo, etcétera). Para el cálculo de las métricas y la discusión de resultados, generamos la siguiente tabla
load("tabla_facil_aprox1_2017-12-10.RData")
head(tabla_facil_aprox1)
En esta tabla se recogen todas las ejecuciones y se ha calculado el valor mínimo y máximo de la función de fitness en cada iteración así como el valor medio. Cada valor está relacionado mediante id.comb (identificador de la combinación) e id.replica (identificador de las pruebas para cada combinación).
Para el cálculo de otras métricas, se genera una tabla en la que se recoge la iteración en la que se haya alcanzado por primera vez el mínimo valor de fitness para cada id.comb e id.replica.
min_iter_facil <- tabla_facil_aprox1 %>%
group_by(id.comb, id.replica) %>%
top_n(-1, min) %>%
top_n(-1, iter)
min_iter_facil
Discusión de resultados
En la siguiente figura se muestra el progreso del algoritmo. Cada una de las líneas representa una ejecución de cada una de las combinaciones de los parámteros.
tabla_facil_aprox1 %>%
ggplot(aes(x = iter,
y = min)) +
geom_line(aes(group = paste(id.comb, id.replica))) +
ylim(0, NA) +
geom_smooth( color = "firebrick")
La curva roja representa un ajuste global mediante GAM (falta referencia) que permite ver la tendencia de la curva. Como suele ser habitual en los algoritmos genéticos, en las primeras iteraciones el algoritmo presenta una fuerte disminución de la función de fitness para después empezar a estabilizarse. En estas $50$ pruebas con $5$ combinaciones de parámetros y $5000$ iteraciones, ninguna ejecución ha alcanzado el óptimo, es decir, una valor de fitness de 0. Al haberse aplicado elitismo en el algoritmo, se aprecian escalones en las líneas. Cabe destacar las línera superiores que se mantienen casi constantes; se corresponden con las ejecuciones realizadas con el valor de probabilidad de mutación igual a 1.
En el siguiente gráfico se puede ver la distribución de los valores de fitness de amnera global
min_iter_facil %>%
ggplot(aes(x = min)) +
geom_histogram()
min_iter_facil %>%
ggplot(aes(x = 'Instancia sencilla' , y = min)) +
geom_jitter(height = 0, width = 0.25, size = 4, alpha = 0.6, color = "steelblue") +
geom_boxplot(alpha = 0)
Además, se muestra que el mínimo se alcanza en r min(min_iter_facil$min), siendo su media y desviación típica, r mean(min_iter_facil$min) y r sd(min_iter_facil$min), respectivamente.
summary(min_iter_facil$min)
sd(min_iter_facil$min)
Ya que nuestro interés se centra en la probabilidad de mutación, podemos comprobar cómo se comporta el aprendizaje del algoritmo dibujando una curva GAM para cada probabilidad de mutación.
tabla_facil_aprox1 %>%
ggplot(aes(x = iter,
y = min,
group = factor(prob_mutacion),
color = prob_mutacion)) +
ylim(0, NA) +
geom_smooth()
En este caso, cuanto más alta es la probabilidad de mutación, el algoritmo tarda más en aprender.
Para enfrentar directamente la probabilidad de mutación y el valor de fitness, cada punto representa el mínimo valor de fitness alcanzado en cada ejecución del algoritmo. Se han agitado los puntos horizontalmente para poder intuir su distribución y se ha dibujado sobre ellos un diagrama de cajas.
min_iter_facil %>%
ggplot(aes(x = factor(prob_mutacion),
y = min)) +
geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) +
geom_boxplot(alpha = 0)
Se aprecia una relación creciente, de hecho la correlación lineal de Pearson obtiene un valor de r cor(min_iter_facil$prob_mutacion, min_iter_facil$min). En el gráfico, no obstante, se aprecia una cierta relación no lineal, con valores de fintess prácticamente indistinguibles para las probabilidades de mutación $0.005$, $0.01$ y $0.05$. Siendo conscientes de una cierta tendencia no lineal, aplicamos aún así un modelo de regresión lineal simple para relacionar la probabilidad de mutación con el fitness obtenido.
summary(lm(formula = min ~ I(prob_mutacion*100), data = min_iter_facil))
En el ajuste lineal, la variable prob_mutacion obtiene un p-valor por debajo de $2e-16$ indicando que hay evidencia estadística significativa de que la probabilidad de mutación afecta al fitness del problema. El parámetro de la regresión para esta variable es de $1.23$, es decir, suponiendo una tendencia lineal, cada incremento del $1%$ en la probabilidad de mutación supondría un incremento de $1.23$ en el valor de fitness.
Por supuesto, hay que destacar el comportamiento anómalo del algoritmo cuando no se considera mutación (prob_mutacion = 0). En este caso los valores son muy superiores al resto de grupos, obteniendo
summary(filter(min_iter_facil, prob_mutacion == 0)$min)
min_iter_facil %>%
ggplot(aes(x = prob_mutacion,
y = iter, size =-min)) +
geom_jitter(height = 0, width = 0.01, color = "steelblue", alpha = 0.7)
min_iter_facil %>%
ggplot(aes(x = factor(prob_mutacion),
y = iter)) +
geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) +
geom_boxplot(alpha = 0)
Se intuye una cierta relación entre la probabilidad de mutación y el número de iteraciones necesarias para alcanzar el mínimo valor del fitness. Valores menores de probabilidades de mutación consiguen encontrar antes, en media, el mínimo valor de la función de fitness.
min_iter_facil %>%
group_by(prob_mutacion) %>%
summarise(mean.min = mean(min),
sd.min = sd(min))
Aunque la media del fitness mínimos para cada probabilidad de mutación fluctúa, la desviación estándar es similar entre los distintos grupos.
¡CONTRASTE DE MEDIAS!
Aproximación 2
En la segunda aproximación, el grid de combinaciones con el que se han realizado las ejecuciones del algoritmo ha sido:
load("sudoko_facil_2017-12-07.RData")
tabla_facil_aprox2 <- analisis_result_tabla(resultados_facil)
saveRDS(tabla_facil_aprox2, "tabla_facil_aprox2.RDS")
tabla_facil_aprox2 <- readRDS("tabla_facil_aprox2.RDS")
min_iter_facil_aprox2 <- tabla_facil_aprox2 %>%
group_by(id.comb, id.replica) %>%
top_n(-1, min) %>%
top_n(-1, iter)
tabla_facil_aprox2 %>%
ggplot(aes(x = iter,
y = min)) +
geom_line(aes(group = paste(id.comb, id.replica))) +
ylim(0, NA) +
geom_smooth( color = "firebrick")
min_iter_facil_aprox2 %>%
filter(min < 60) %>%
ggplot(aes(x = prob_mutacion,
y = min)) +
geom_point(size = 3, color = "steelblue", alpha = 0.7) +
geom_smooth()
min_iter_facil_aprox2 %>%
ggplot(aes(x = prob_mutacion,
y = iter)) +
geom_point( size = 3, color = "steelblue", alpha = 0.7)
Instancia compleja
Aproximación 1
La discusión de resultados seguirá la misma línea que lo expuesto en la sección BLA por lo que se omitirá el detalle de los cálculos.
Para la instancia compleja, se ha elegido un tamaño de población de $100$.
Los resultados de esta ejecución se encuentran en el archivo actividad/tabla_dificil_aprox1_2017-12-11.RData
load("tabla_dificil_aprox1_2017-12-11.RData")
min_iter_dificil <- tabla_dificil_aprox1 %>%
group_by(id.comb, id.replica) %>%
top_n(-1, min) %>%
top_n(-1, iter)
En la siguiente figura se muestra el progreso del algoritmo. Cada una de las líneas representa una ejecución de cada una de las combinaciones de los parámetros.
tabla_dificil_aprox1 %>%
ggplot(aes(x = iter,
y = min)) +
geom_line(aes(group = paste(id.comb, id.replica))) +
ylim(0, NA) +
geom_smooth( color = "firebrick")
La tendencia es similar a la seguida en la instancia sencilla, con una bajada de la función de fitness muy acusada en las primeras iteraciones y una estabilización general después.
Igualmente parece existir una relación entre la probabilidad de mutación y la convergencia del algoritmo de forma que con un nivel menos de probabilidad de mutación se obtiene una convergencia más rápida y, con el líminte impuesto de $5000$ iteraciones, llega a un menor valor de la función de fitness las probabilidades de mutación menores.
tabla_dificil_aprox1 %>%
ggplot(aes(x = iter,
y = min,
group = factor(prob_mutacion),
color = prob_mutacion)) +
ylim(0, NA) +
geom_smooth()
Para enfrentar directamente la probabilidad de mutación y el valor de fitness, cada punto representa el mínimo valor de fitness alcanzado en cada ejecución del algoritmo. Se han agitado los puntos horizontalmente para poder intuir su distribución y se ha dibujado sobre ellos un diagrama de cajas.
min_iter_dificil <- tabla_dificil_aprox1 %>%
group_by(id.comb, id.replica) %>%
top_n(-1, min) %>%
top_n(-1, iter)
min_iter_dificil %>%
ggplot(aes(x = prob_mutacion,
y = min)) +
geom_point(size = 3, color = "steelblue", alpha = 0.7) +
geom_smooth()
min_iter_dificil %>%
ggplot(aes(x = factor(prob_mutacion),
y = min)) +
geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) +
geom_boxplot(alpha = 0)
De nuevo se aprecia una relación, esta vez más claramente no lineal. En la instancia sencilla, el resultado medio de la función de fitness para las probabilidades de mutación de $0.01$ y $0.05$ era prácticamente indistinguibles. Sin embargo, en la instancia compleja el mínimo de la función de fitness (en media) se alcanza con una probabilidad de mutación de $0.05$.
Se aprecia una relación creciente. La correlación lineal de Pearson obtiene un valor de r cor(min_iter_dificil$prob_mutacion, min_iter_dificil$min). Y, si aplicamos un modelo de regresión lineal simple
summary(lm(formula = min ~ I(prob_mutacion*100), data = min_iter_dificil))
En el ajuste lineal, la variable prob_mutacion obtiene un p-valor por debajo de $2e-16$ indicando que hay evidencia estadística significativa de que la probabilidad de mutación afecta al fitness del problema. El parámetro de la regresión para esta variable es de $1.18597$, es decir, suponiendo una tendencia lineal, cada incremento del $1%$ en la probabilidad de mutación supondría un incremento de $1.18597$ en el valor de fitness.
min_iter_dificil %>%
ggplot(aes(x = factor(prob_mutacion),
y = iter)) +
geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) +
geom_boxplot(alpha = 0)
Aproximación 2
load("sudoko_dificil_2017-12-11.RData")
tabla_dificil_aprox2 <- analisis_result_tabla(resultados_dificil)
saveRDS(tabla_dificil_aprox2, "tabla_dificil_aprox2.RDS")
tabla_dificil_aprox2 <- readRDS("tabla_dificil_aprox2.RDS")
min_iter_dificil_aprox2 <- tabla_dificil_aprox2 %>%
group_by(id.comb, id.replica) %>%
top_n(-1, min) %>%
top_n(-1, iter)
tabla_dificil_aprox2 %>%
ggplot(aes(x = iter,
y = min)) +
geom_line(aes(group = paste(id.comb, id.replica))) +
ylim(0, NA) +
geom_smooth( color = "firebrick")
min_iter_dificil_aprox2 %>%
filter(min < 60) %>%
ggplot(aes(x = prob_mutacion,
y = min)) +
geom_point(size = 3, color = "steelblue", alpha = 0.7) +
geom_smooth()
min_iter_dificil_aprox2 %>%
ggplot(aes(x = prob_mutacion,
y = iter)) +
geom_point( size = 3, color = "steelblue", alpha = 0.7)
Comparación
min_iter_facil$instancia <- 'sencilla'
min_iter_dificil$instancia <- 'compleja'
min_iter <- bind_rows(min_iter_facil, min_iter_dificil)
min_iter %>%
filter(prob_mutacion != 0) %>%
ggplot(aes(x = instancia, y = min)) +
geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) +
geom_boxplot(alpha = 0)
min_iter %>%
filter(prob_mutacion != 0) %>%
ggplot(aes(x = instancia, y = iter)) +
geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) +
geom_boxplot(alpha = 0)
min_iter %>%
filter(prob_mutacion != 0) %>%
ggplot(aes(x = prob_mutacion, y = min, color = instancia)) +
geom_point(size = 3, alpha = 0.7)
min_iter %>%
filter(prob_mutacion != 0) %>%
ggplot(aes(x = prob_mutacion, y = min, color = instancia)) +
geom_smooth(size = 3, alpha = 0.7)
min_iter %>%
filter(prob_mutacion != 0) %>%
ggplot(aes(x = prob_mutacion, y = iter, color = instancia)) +
geom_smooth(size = 3, alpha = 0.7)
El objetivo último es el de evaluar cómo impactan los diferentes valores de la probabilidad de mutación en el algoritmo genético.
Para ello se ha seguido la siguiente metodología.
- Probar el algoritmo con dos instancias: una compleja y otra sencilla.
- Con una configuración de parámetros fija, variar la probabilidad de mutación para ver cómo ésta puede afectar al desarrollo del algoritmo.
- Construir 50 muestras aleatorias de los parámetros para poder cubrir el mayor espacio y poder extraer conclusiones del algoritmo.
Los parámetros que rigen el algoritmo genético que se van a variar son:
- Tamaño de la población.
- En la selección de padres por torneo, el tamaño del torneo.
- Probabilidad de cruce
- Probabilidad de mutación.
Para evaluar la idoneidad del algoritmo genético, se han generado $100$ combinaciones aleatorias de los parámetros del algoritmo.
- Tamaño de la población.: entre $10$ y $200$ individuos.
- Tamaño del torneo: entre $2$ y $10$.
- Probabilidad de cruce: valor con distrbución uniforme $U(0.5, 1)$.
- Probabilidad de mutación: valor con distrbución uniforme $U(0.01, 0.2)$.
library(dplyr)
set.seed(1234)
n <- 100
tam_poblacion <- sample(10:200, size = n, replace = TRUE)
prob_cruce <- runif(n = n, min = 0.5, max = 1)
tam_torneo <- sample(2:10, size = n, replace =TRUE)
prob_mutacion <- runif(n = n, min = 0.001, max = 0.2)
grid <- data_frame(tam_poblacion,
prob_cruce,
tam_torneo,
prob_mutacion)
head(grid)
Cada combinación de parámetros se ha replicado $10$ veces para poder obtener datos sobre los que realizar análisis estadísticos.
La metodología que se ha llevado ha sido la de lanzar
Conclusiones
La resolución del Sudoku mediante el algoritmo genético implementado no parece ser la técnica más adecuada. En muchas de la mayoría de las combinaciones, el problema del Sudoku no estaba resuelto en las $5000$ primeras iteraciones, por lo que no se ha observado una eficiencia.
Una parte que se podría tratar de estudiar es la de si la representación de los individuos mediante un gen en cada celda con un número comprendido entre 1 y 9 es la más adecuada. Se propone como trabajo futuro el estudio de otras respresentaciones, como la representación con permutaciones. Es decir, cada individuo estaría representado por un genotipo
$$
\begin{aligned}
s = (&1, 2, \ldots, 9,\
&1, 2, \ldots, 9,\
&1, 2, \ldots, 9,\
&\ldots,\
&1, 2, \ldots, 9)
\end{aligned}
$$
Y cualquier individuo debería ser una permutación, lo que supondría elegir otro tipo de mutación que la random resetting elegida para esta implementación.
Referencias
papabloblo/sudoku documentation built on May 30, 2019, 3:45 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(sudoku) library(tidyverse)
Introducción
Los algoritmos propios de la computación evolutiva pueden aplicarse en campos muy diversos. El objetivo de esta actividad es el de describir el comportamiento de los algoritmos genéticos aplicados a la resolución del pasatiempo conocido como Sudoku.
Comenzaremos describiendo en qué consiste el problema del sudoku y definiéndolo de una manera formal. Después, se comentarácómo se ha llevado a cabo la implementación en R para pasar a hacer una evaluación experimental que nos permita discutir el comportamiento.
Descripción del problema del sudoku
El juego del Sudoku comenzó a ser popular en Japón en 1986 y consiguió popularidad internacional en 2005 [@Lynce2006]. Se trata de un pasatiempo cuyas reglas son simples de entender pero su resolución puede ser todo un reto. Además, un sudoku bien construido tiene una única solución y puede ser resuelto mediante el razonamiento sin tener que recurrir a procedimientos de ensayo y error.
Los sudokus suelen tener asociado un nivel de dificultad que suele depender del número de casillas rellenas en el inicio. En [@McGuire2014] se muestra que el número mínimo de casillas que deben de estar rellenas para que un sudoku tenga solución (única) es 17.
¿En qué consiste?
El juego del Sudoku tradicional se representa por una matriz de dimensión $9x9$ en la que se disponen números enteros del 1 al 9. En un inicio solo algunas de las casillas están ocupadas con números siendo el objetivo del juego rellenar las casillas vacías.
Los elementos principales del sudoku serán las filas, las columnas y las subcuadrículas como se muestra en (imagen).
Para que el sudoku se considere resuelto, se ha de verificar que:
- cada fila contiene una permutación de los 9 enteros apareciendo cada uno de ellos una y solo una vez,
- cada columna contiene una permutación de los 9 enteros apareciendo cada uno de ellos una y solo una vez,
- cada subcuadrícula debe contener una permutación de los 9 enteros apareciendo cada uno de ellos una y solo una vez.
Estas reglas de resolución explican porqué Sudoku significa "un solo número" en japonés [@Lynce2006].
Formalmente
El problema puede ser expresado matemáticamente de múltiples formas. El desarrollo elegido tiene la motivación de que la aplicación del algoritmo evolutivo sea más natural.
El sudoku se se expresa como un vector $\mathbf{s} = (s_1, s_2, \ldots, s_{81})$ donde $s_i \in {0, 1, 2,\ldots,9}$ para $i = 1, \ldots, 81$. Los elementos sobre los que se deben evaluar la correcta resolución del problema y que servirán para calcular la función fitness del algoritmo evolutivo se describen matemáticamente como
$$ columna(j) = {i \in {1,\ldots, 81}\, | \, mod(i, 9)= mod(j, 9)} $$
$$ fila(j) = {i \in {1,\ldots, 81}\, | \, \lfloor j/9 \rfloor = \lfloor i/9\rfloor} $$
donde $\lfloor a \rfloor$ se refiere a la parte entera de $a$.
Por otro lado, definimos la funciones:
$$ c(i) = \sum_{{j \in columna(i)|j \ne i}} I(s_i = s_j) $$ $$ f(i) = \sum_{{j \in fila(i)|j \ne i}} I(s_i = s_j) \ $$ $$ s(i) = \sum_{{j \in subcuadrícula(i)|j \ne i}} I(s_i = s_j) \ $$
El sudoku se considerará resuelto si y solo si $\sum_i f(i) = 0$, $\sum_i c(i) = 0$ y $\sum_i s(i) = 0$, o lo que es lo mismo
$$ \sum_{i=1}^{81} \left[f(i) + c(i) + s(i)\right] = 0 $$ Podemos definir la siguiente función
$$ fitness(s) = \sum_{i=1}^{81} f(i) + c(i) + s(i) $$
Cualquier configuración de números que no resuelva el sudoku tendrá un $$fitness(s) > 0$$ por lo que, la solución que resuelve el sudoku se puede expresar como $$ s^* = \underset{s}{\arg\min} \, fitness(s) $$
-->
Implementación
El juego del sudoku se va a abordar mediante un algoritmo genético. Hay distintas estrategias que se podrían seguir para este problema. En este caso, la implementación se ha realizado conforme a las características que se detallan a continuación.
Estructura de códigos en R
Los códigos implementados se encuentran en la carpeta \R.
Para que los resultados de las siguientes secciones sea reproducibles, se fija la semilla para la generación de valores aleatorios
set.seed(123)
Representación de cada individuo
Cada celda $s_i$ del sudoku estará asociada a un gen representado por números enteros del 1 al 9, es decir, $s_i\in {1, 2, \ldots, 9}$. El genotipo estará entonces representado por una cadena de enteros. Lo que es lo mismo, se realizará una muestra con reemplazamiento de los enteros del 1 al 9 con una longitud determinada.
La función random_integer_representation() es la que permite obtener este tipo de representación. Necesita los argumentos:
valores: vector con los valores posibles que puede tomar cada gen.tam: número entero con el tamaño del genotipo (número de genes).
Como ejemplo:
random_integer_representation(valores = 1:9, tam = 9)
Aunque el genotipo de cada individuo contendrá un total de $81$ genes, solo se aplicará el algoritmo sobre aquellos genes (casillas) que estén vacíos en la solución inicial. Por tanto, el algoritmo genético manejará el genotipo $s_0$ de forma que
$$ s_i = s^1_i si $$
Función de adaptación (fitness)
La función de fitness que el algoritmo tratará de minimizar será la que aparece en la ecuación ]bla]. Como se dijo anteriormente, un valor de $fitness(s)=0$ asegura que el juego ha sido resuelto. El valor máximo de esta función se alcanzará cuando todas las casillas (genes) contengan el mismo valor y tendrá un valor de $fitness(s) = 972$
Las funciones implementadas para el cálculo de la función de fitness se encuentran en \R\fitness.R.
La función fitness_sudoku() calcula la función de fitness para un individuo dado y necesita los argumentos:
x: un vector con el genotipo del individuo.ind_cuadricula: un vector indicando qué posiciones equivalen a la representación de las subcuadrículas.solucion_inicial: un vector representando la solución inicial del sudoku que contendráNAen aquellas posiciones que originalmente aparecían vacías.
En nuestro caso, el ind_cuadricula se representa como
ind_cuadricula <- c( rep(rep(x = 1:3, each = 3), 3), rep(rep(x = 4:6, each = 3), 3), rep(rep(x = 7:9, each = 3), 3)) ind_cuadricula
Viendo este vector como una matriz, se pone de relieve cómo están representadas las subcuadrículas
matrix(ind_cuadricula, nrow = 9, ncol = 9)
Un ejemplo de solución inicial del sudoku se podría representar así:
sudoku_facil <- c(NA, NA, NA, NA, 3, NA, NA, NA, 4, NA, 9, NA, 4, NA, 6, NA, 7, NA, NA, 5, NA, NA, NA, NA, 3, 8, NA, NA, NA, NA, NA, 7, 8, NA, NA, 3, 3, NA, NA, NA, NA, NA, 6, 9, NA, 5, 4, NA, 6, NA, NA, NA, 2, NA, 7, NA, 5, NA, 2, 4, NA, NA, NA, 9, 8, 4, NA, 6, 5, 2, NA, NA, NA, 2, 6, NA, 8, NA, NA, NA, 9)
El genotipo de los individuos estará representado por el número de casillas vacías en la solución inicial, es decir,
tam_individuos <- sum(is.na(sudoku_facil)) tam_individuos
Simulamos un individuo con 48 genes, en este caso,
individuo <- random_integer_representation(valores = 1:9, tam = tam_individuos)
Ya podemos calcular el valor de la función de fitness para este individuo.
fitness_sudoku(x = individuo, ind_cuadricula = ind_cuadricula, solucion_inicial = sudoku_facil)
Inicialización de la población
La población inicial del algoritmo se realizatá de forma aleatoria, es decir, cada individuo estará representado por los genes que le corresponda y cada gen se inicializará como un número entero aleatorio entre 1 y 9.
Esto se permite gracias a la función generacion_poblacion() cuyos argumentos son:
valores_posibles: vector con los valores posibles que puede tomar cada gen.num_genes: número entero con el tamaño del genotipo (número de genes).tam_poblacion: entero con el número de individuos de los que estará formada la población.
El resultado es una lista compuesta por los distintos individuos. Por ejemplo,
poblacion <- generacion_poblacion(valores_posibles = 1:9, num_genes = tam_individuos, tam_poblacion = 3) poblacion
Selección de padres
Para poder realizar el cruce entre individuos que se describe en la siguiente sección, es necesario elegir los padres, es decir, aquellos individuos entre los que realizar el cruce. En este caso se realiza la selección por torneo.
Tenemos la función seleccion_padres() cuyos argumentos son:
num_padres:poblacion: lista con los individuos que compoennfitness_poblacion: vector con el valor de fitness para cada genotipo de la población.k: entero con el número de competidores del torneo
El resultado es una lista con los padres elegidos por torneo.
padres <- seleccion_padres(num_padres = 2, poblacion = poblacion, fitness_poblacion = sapply(poblacion, fitness_sudoku, ind_cuadricula = ind_cuadricula, solucion_inicial = sudoku_facil), k = 2) padres
Cruce
El cruce entre los padres se realizará mediante cruce por un punto.
Función one_point_crossover() con argumentos:
padres: lista con los padres seleccionados para el cruce.prob_cruce:numericen [0, 1] representando la probabilidad de que los padres se crucen y tengan descendencia o se mantengan.
resultado_cruce <- one_point_crossover(padres = padres, prob_cruce = 0.75) resultado_cruce
Mutación
Una vez producido el cruce entre padres y obtenida su descendencia, se aplicará mutación random resetting, es decir, cada uno de los genes de cada individuo podrá sufrir una mutación con una probabilidad dada. En el caso de que se produzca una mutación en un gen, esta podrá tomar un valor entero aleatorio entre 1 y 9.
random_resetting() con parámetros:
x: individuo sobre el que provocar mutacionesprob: probabilidad de que se aplique mutación en cada gen.valores_posibles: vector con valores posibles para la mutación de cada gen.
random_resetting(individuo, prob = 0.1, valores_posibles = 1:9)
Selección de supervivientes
Se utilizará un modelo de tipo generacional, reemplazando la población actual enteramente por una nueva aplicando, eso sí, elitismo para conservar el emjor individuo en la siguiente generación y eliminando el peor.
El algoritmo
Una ejecución del algoritmo se realizará gracias a la función algoritmo_genetico() con parámetros
poblacion_inicial: lista con los individuos que compone la población inicial.funcion_fitness: función que calculael fitness de cada individuo de la población.genes_fijos: ...valores_mutacion: vector con los valores posibles para la mutación.num_padres: entero con el número de padres.prob_cruce: probabilidad de que se produzca cruce entre padres.tam_torneo: tamaño del torneo en la selección de padres por torneo.prob_mutacion: probabilidad de que se produzca una mutación en cada gen del individuo.tam_poblacion: tamaño de la población ¿?max_iter: iteraciones máximas que realizará el algoritmoprint_each: imprimir por pantalla información del algoritmo cadaprint_eachiteraciones.
Se ha elegido como lenguaje de programación R. R es un tipo de lenguaje Open Source con una gran aceptación dentro del área del data science.
Todas las funciones se han implementado de tal forma que puedan utilizarse para aplicar un algoritmo genético en general y no se han particularizado para el problema del Sudoku a excepción de la función de fitness.
Toda la documentación de cada una de las funciones implementadas se encuentra en el pdf adjunto.
Evaluación experimental
Los algoritmos genéticos llevan asociados una serie de parámetros que rigen su comportamiento y pueden impactar en los resultados obtenidos. Por lo tanto, es fundamental conocer cómo afectan los distintos parámetros al algoritmo para poder tener un conocimiento sobre los mismos que nos permita ser conscientes de los cambios. Por ejemplo, en función del problema, nos puede interesar favorecer una combinación de parámetros que necesite un menor tiempo computacional (un menor número de iteraciones) para obtener una solución cuya función de fitness sea ligeramente peor (en media) que la de otra combinación de parámetros pero que necesite un alto número de iteraciones.
Metodología
Para contrastar el comportamiento del algoritmo genético aplicado al problema de resolver un Sudoku, se han elegido dos configuraciones iniciales, una fácil con 33 casillas iniciales rellenas y otra difícil, con 30 casillas rellenas. Como se dijo en la sección BLA, se considera que la dificultad del problema del Sudoku es inversamente proporcional al número de celdas rellenas al principio.
En la implementación descrita en la sección BLA del algoritmo genético, los parámetros del algoritmo serían:
- Tamaño de la población.
- En la selección de padres por torneo, el tamaño del torneo.
- Probabilidad de cruce
- Probabilidad de mutación.
En nuestro caso, una de las preguntas fundamentales que queremos responder es cómo impactan los diferentes valores de la probabilidad de mutación en el algoritmo genético.
A continuación se describen las dos metodologías seguidas.
Aproximación 1
Se han fijado todos los parámetros del algoritmo a excepción de la probabilidad de mutación para la cual se han elegido valores. Se ha ejecutado el algoritmo 50 veces.
Aproximación 2
Para tratar de obtener un resultado del impacto global de los parámetros, se han generado $50$ configuraciones aleatorias de los parámetros con la siguiente lógica: - El tamaño de la población se fija en $10$ para la instancia sencilla y $100$ para la instancia compleja. - Probabilidad de cruce se supone una variable aleatoria uniforme $U(0.5, 1)$. - Tamaño del torneo: valor entero aleatorio entre $2$ y $10$, ambos incluidos. - Probabilidad de cruce: valor con distribución uniforme $U(0.5, 1)$. - Probabilidad de mutación: permutación aleatoria de 50 números igualmente espaciados entre 0 y 0.25. - $10$ ejecuciones para cada combinación de parámetros
Métricas y gráficas
Instancia sencilla
La instancia sencilla se ha tomado de http://www.sudoku-online.org/ en la categoría fácil 10 (#1217). Su vector se implementa de la siguiente forma
sudoku_facil <- c(NA, NA, NA, NA, 3, NA, NA, NA, 4, NA, 9, NA, 4, NA, 6, NA, 7, NA, NA, 5, NA, NA, NA, NA, 3, 8, NA, NA, NA, NA, NA, 7, 8, NA, NA, 3, 3, NA, NA, NA, NA, NA, 6, 9, NA, 5, 4, NA, 6, NA, NA, NA, 2, NA, 7, NA, 5, NA, 2, 4, NA, NA, NA, 9, 8, 4, NA, 6, 5, 2, NA, NA, NA, 2, 6, NA, 8, NA, NA, NA, 9)
Nota: tanto la instancia compleja, como la instancia sencilla se recogen el script \sudoku_prueba.R en los objetos sudoku_facil y sudoku_dificil, respectivamente.
Aproximación 1
Se genera un grid con la siguiente configuración de parámetros
grid <- data.frame(tam_poblacion = 10, prob_cruce = 0.9, tam_torneo = 2, prob_mutacion = c(0, 0.01, 0.05, 0.1, 0.15, 0.25)) grid
Como ya se comentó, para esta primera aproximación se mantienen todos los parámetros constantes a excepción de la probabilidad de mutación que toma valores $0$, $0.01$, $0.05$, $0.1$, $0.15$ y $0.25$.
Para cada combinación de parámetros se ha repetido la ejecución $50$ veces dejando un máximo de $5000$ iteraciones.
Esta ejecución se recoge en
library(parallel) set.seed(1234) ini <- Sys.time() resultados_dificil <- mcmapply(pruebas_ga, num_pruebas = 50, num_padres = 2, tam_poblacion = grid$tam_poblacion, prob_cruce = grid$prob_cruce, tam_torneo = grid$tam_torneo, prob_mutacion = grid$prob_mutacion, max_iter = 5000, print_each = 250, MoreArgs = list(generacion_poblacion_ini = generacion_poblacion, genes_fijos = sudoku_facil, funcion_fitness = fitness_sudoku, valores_posibles = 1:9, valores_mutacion = 1:9 ) , SIMPLIFY = FALSE ) fin <- Sys.time() fin-ini save(resultados_dificil, file = paste0("data/prueba_facil_", Sys.Date(), ".RData"))
Los resultados de esta ejecución se encuentran en el archivo actividad/prueba_facil_2017-12-10.RData. El objeto resultados_facil contiene una gran cantidad de información (población de cada iteración, fitness de cada individuo, etcétera). Para el cálculo de las métricas y la discusión de resultados, generamos la siguiente tabla
load("tabla_facil_aprox1_2017-12-10.RData") head(tabla_facil_aprox1)
En esta tabla se recogen todas las ejecuciones y se ha calculado el valor mínimo y máximo de la función de fitness en cada iteración así como el valor medio. Cada valor está relacionado mediante id.comb (identificador de la combinación) e id.replica (identificador de las pruebas para cada combinación).
Para el cálculo de otras métricas, se genera una tabla en la que se recoge la iteración en la que se haya alcanzado por primera vez el mínimo valor de fitness para cada id.comb e id.replica.
min_iter_facil <- tabla_facil_aprox1 %>% group_by(id.comb, id.replica) %>% top_n(-1, min) %>% top_n(-1, iter) min_iter_facil
Discusión de resultados
En la siguiente figura se muestra el progreso del algoritmo. Cada una de las líneas representa una ejecución de cada una de las combinaciones de los parámteros.
tabla_facil_aprox1 %>% ggplot(aes(x = iter, y = min)) + geom_line(aes(group = paste(id.comb, id.replica))) + ylim(0, NA) + geom_smooth( color = "firebrick")
La curva roja representa un ajuste global mediante GAM (falta referencia) que permite ver la tendencia de la curva. Como suele ser habitual en los algoritmos genéticos, en las primeras iteraciones el algoritmo presenta una fuerte disminución de la función de fitness para después empezar a estabilizarse. En estas $50$ pruebas con $5$ combinaciones de parámetros y $5000$ iteraciones, ninguna ejecución ha alcanzado el óptimo, es decir, una valor de fitness de 0. Al haberse aplicado elitismo en el algoritmo, se aprecian escalones en las líneas. Cabe destacar las línera superiores que se mantienen casi constantes; se corresponden con las ejecuciones realizadas con el valor de probabilidad de mutación igual a 1.
En el siguiente gráfico se puede ver la distribución de los valores de fitness de amnera global
min_iter_facil %>% ggplot(aes(x = min)) + geom_histogram()
min_iter_facil %>% ggplot(aes(x = 'Instancia sencilla' , y = min)) + geom_jitter(height = 0, width = 0.25, size = 4, alpha = 0.6, color = "steelblue") + geom_boxplot(alpha = 0)
Además, se muestra que el mínimo se alcanza en r min(min_iter_facil$min), siendo su media y desviación típica, r mean(min_iter_facil$min) y r sd(min_iter_facil$min), respectivamente.
summary(min_iter_facil$min) sd(min_iter_facil$min)
Ya que nuestro interés se centra en la probabilidad de mutación, podemos comprobar cómo se comporta el aprendizaje del algoritmo dibujando una curva GAM para cada probabilidad de mutación.
tabla_facil_aprox1 %>% ggplot(aes(x = iter, y = min, group = factor(prob_mutacion), color = prob_mutacion)) + ylim(0, NA) + geom_smooth()
En este caso, cuanto más alta es la probabilidad de mutación, el algoritmo tarda más en aprender.
Para enfrentar directamente la probabilidad de mutación y el valor de fitness, cada punto representa el mínimo valor de fitness alcanzado en cada ejecución del algoritmo. Se han agitado los puntos horizontalmente para poder intuir su distribución y se ha dibujado sobre ellos un diagrama de cajas.
min_iter_facil %>% ggplot(aes(x = factor(prob_mutacion), y = min)) + geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) + geom_boxplot(alpha = 0)
Se aprecia una relación creciente, de hecho la correlación lineal de Pearson obtiene un valor de r cor(min_iter_facil$prob_mutacion, min_iter_facil$min). En el gráfico, no obstante, se aprecia una cierta relación no lineal, con valores de fintess prácticamente indistinguibles para las probabilidades de mutación $0.005$, $0.01$ y $0.05$. Siendo conscientes de una cierta tendencia no lineal, aplicamos aún así un modelo de regresión lineal simple para relacionar la probabilidad de mutación con el fitness obtenido.
summary(lm(formula = min ~ I(prob_mutacion*100), data = min_iter_facil))
En el ajuste lineal, la variable prob_mutacion obtiene un p-valor por debajo de $2e-16$ indicando que hay evidencia estadística significativa de que la probabilidad de mutación afecta al fitness del problema. El parámetro de la regresión para esta variable es de $1.23$, es decir, suponiendo una tendencia lineal, cada incremento del $1%$ en la probabilidad de mutación supondría un incremento de $1.23$ en el valor de fitness.
Por supuesto, hay que destacar el comportamiento anómalo del algoritmo cuando no se considera mutación (prob_mutacion = 0). En este caso los valores son muy superiores al resto de grupos, obteniendo
summary(filter(min_iter_facil, prob_mutacion == 0)$min)
min_iter_facil %>% ggplot(aes(x = prob_mutacion, y = iter, size =-min)) + geom_jitter(height = 0, width = 0.01, color = "steelblue", alpha = 0.7)
min_iter_facil %>% ggplot(aes(x = factor(prob_mutacion), y = iter)) + geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) + geom_boxplot(alpha = 0)
Se intuye una cierta relación entre la probabilidad de mutación y el número de iteraciones necesarias para alcanzar el mínimo valor del fitness. Valores menores de probabilidades de mutación consiguen encontrar antes, en media, el mínimo valor de la función de fitness.
min_iter_facil %>% group_by(prob_mutacion) %>% summarise(mean.min = mean(min), sd.min = sd(min))
Aunque la media del fitness mínimos para cada probabilidad de mutación fluctúa, la desviación estándar es similar entre los distintos grupos.
¡CONTRASTE DE MEDIAS!
Aproximación 2
En la segunda aproximación, el grid de combinaciones con el que se han realizado las ejecuciones del algoritmo ha sido:
load("sudoko_facil_2017-12-07.RData") tabla_facil_aprox2 <- analisis_result_tabla(resultados_facil) saveRDS(tabla_facil_aprox2, "tabla_facil_aprox2.RDS")
tabla_facil_aprox2 <- readRDS("tabla_facil_aprox2.RDS") min_iter_facil_aprox2 <- tabla_facil_aprox2 %>% group_by(id.comb, id.replica) %>% top_n(-1, min) %>% top_n(-1, iter)
tabla_facil_aprox2 %>% ggplot(aes(x = iter, y = min)) + geom_line(aes(group = paste(id.comb, id.replica))) + ylim(0, NA) + geom_smooth( color = "firebrick")
min_iter_facil_aprox2 %>% filter(min < 60) %>% ggplot(aes(x = prob_mutacion, y = min)) + geom_point(size = 3, color = "steelblue", alpha = 0.7) + geom_smooth()
min_iter_facil_aprox2 %>% ggplot(aes(x = prob_mutacion, y = iter)) + geom_point( size = 3, color = "steelblue", alpha = 0.7)
Instancia compleja
Aproximación 1
La discusión de resultados seguirá la misma línea que lo expuesto en la sección BLA por lo que se omitirá el detalle de los cálculos.
Para la instancia compleja, se ha elegido un tamaño de población de $100$.
Los resultados de esta ejecución se encuentran en el archivo actividad/tabla_dificil_aprox1_2017-12-11.RData
load("tabla_dificil_aprox1_2017-12-11.RData")
min_iter_dificil <- tabla_dificil_aprox1 %>% group_by(id.comb, id.replica) %>% top_n(-1, min) %>% top_n(-1, iter)
En la siguiente figura se muestra el progreso del algoritmo. Cada una de las líneas representa una ejecución de cada una de las combinaciones de los parámetros.
tabla_dificil_aprox1 %>% ggplot(aes(x = iter, y = min)) + geom_line(aes(group = paste(id.comb, id.replica))) + ylim(0, NA) + geom_smooth( color = "firebrick")
La tendencia es similar a la seguida en la instancia sencilla, con una bajada de la función de fitness muy acusada en las primeras iteraciones y una estabilización general después.
Igualmente parece existir una relación entre la probabilidad de mutación y la convergencia del algoritmo de forma que con un nivel menos de probabilidad de mutación se obtiene una convergencia más rápida y, con el líminte impuesto de $5000$ iteraciones, llega a un menor valor de la función de fitness las probabilidades de mutación menores.
tabla_dificil_aprox1 %>% ggplot(aes(x = iter, y = min, group = factor(prob_mutacion), color = prob_mutacion)) + ylim(0, NA) + geom_smooth()
Para enfrentar directamente la probabilidad de mutación y el valor de fitness, cada punto representa el mínimo valor de fitness alcanzado en cada ejecución del algoritmo. Se han agitado los puntos horizontalmente para poder intuir su distribución y se ha dibujado sobre ellos un diagrama de cajas.
min_iter_dificil <- tabla_dificil_aprox1 %>% group_by(id.comb, id.replica) %>% top_n(-1, min) %>% top_n(-1, iter)
min_iter_dificil %>% ggplot(aes(x = prob_mutacion, y = min)) + geom_point(size = 3, color = "steelblue", alpha = 0.7) + geom_smooth()
min_iter_dificil %>% ggplot(aes(x = factor(prob_mutacion), y = min)) + geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) + geom_boxplot(alpha = 0)
De nuevo se aprecia una relación, esta vez más claramente no lineal. En la instancia sencilla, el resultado medio de la función de fitness para las probabilidades de mutación de $0.01$ y $0.05$ era prácticamente indistinguibles. Sin embargo, en la instancia compleja el mínimo de la función de fitness (en media) se alcanza con una probabilidad de mutación de $0.05$.
Se aprecia una relación creciente. La correlación lineal de Pearson obtiene un valor de r cor(min_iter_dificil$prob_mutacion, min_iter_dificil$min). Y, si aplicamos un modelo de regresión lineal simple
summary(lm(formula = min ~ I(prob_mutacion*100), data = min_iter_dificil))
En el ajuste lineal, la variable prob_mutacion obtiene un p-valor por debajo de $2e-16$ indicando que hay evidencia estadística significativa de que la probabilidad de mutación afecta al fitness del problema. El parámetro de la regresión para esta variable es de $1.18597$, es decir, suponiendo una tendencia lineal, cada incremento del $1%$ en la probabilidad de mutación supondría un incremento de $1.18597$ en el valor de fitness.
min_iter_dificil %>% ggplot(aes(x = factor(prob_mutacion), y = iter)) + geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) + geom_boxplot(alpha = 0)
Aproximación 2
load("sudoko_dificil_2017-12-11.RData") tabla_dificil_aprox2 <- analisis_result_tabla(resultados_dificil) saveRDS(tabla_dificil_aprox2, "tabla_dificil_aprox2.RDS")
tabla_dificil_aprox2 <- readRDS("tabla_dificil_aprox2.RDS") min_iter_dificil_aprox2 <- tabla_dificil_aprox2 %>% group_by(id.comb, id.replica) %>% top_n(-1, min) %>% top_n(-1, iter)
tabla_dificil_aprox2 %>% ggplot(aes(x = iter, y = min)) + geom_line(aes(group = paste(id.comb, id.replica))) + ylim(0, NA) + geom_smooth( color = "firebrick")
min_iter_dificil_aprox2 %>% filter(min < 60) %>% ggplot(aes(x = prob_mutacion, y = min)) + geom_point(size = 3, color = "steelblue", alpha = 0.7) + geom_smooth()
min_iter_dificil_aprox2 %>% ggplot(aes(x = prob_mutacion, y = iter)) + geom_point( size = 3, color = "steelblue", alpha = 0.7)
Comparación
min_iter_facil$instancia <- 'sencilla' min_iter_dificil$instancia <- 'compleja' min_iter <- bind_rows(min_iter_facil, min_iter_dificil)
min_iter %>% filter(prob_mutacion != 0) %>% ggplot(aes(x = instancia, y = min)) + geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) + geom_boxplot(alpha = 0)
min_iter %>% filter(prob_mutacion != 0) %>% ggplot(aes(x = instancia, y = iter)) + geom_jitter(height = 0, size = 3, color = "steelblue", alpha = 0.7) + geom_boxplot(alpha = 0)
min_iter %>% filter(prob_mutacion != 0) %>% ggplot(aes(x = prob_mutacion, y = min, color = instancia)) + geom_point(size = 3, alpha = 0.7)
min_iter %>% filter(prob_mutacion != 0) %>% ggplot(aes(x = prob_mutacion, y = min, color = instancia)) + geom_smooth(size = 3, alpha = 0.7)
min_iter %>% filter(prob_mutacion != 0) %>% ggplot(aes(x = prob_mutacion, y = iter, color = instancia)) + geom_smooth(size = 3, alpha = 0.7)
El objetivo último es el de evaluar cómo impactan los diferentes valores de la probabilidad de mutación en el algoritmo genético.
Para ello se ha seguido la siguiente metodología.
- Probar el algoritmo con dos instancias: una compleja y otra sencilla.
- Con una configuración de parámetros fija, variar la probabilidad de mutación para ver cómo ésta puede afectar al desarrollo del algoritmo.
- Construir 50 muestras aleatorias de los parámetros para poder cubrir el mayor espacio y poder extraer conclusiones del algoritmo.
Los parámetros que rigen el algoritmo genético que se van a variar son:
- Tamaño de la población.
- En la selección de padres por torneo, el tamaño del torneo.
- Probabilidad de cruce
- Probabilidad de mutación.
Para evaluar la idoneidad del algoritmo genético, se han generado $100$ combinaciones aleatorias de los parámetros del algoritmo.
- Tamaño de la población.: entre $10$ y $200$ individuos.
- Tamaño del torneo: entre $2$ y $10$.
- Probabilidad de cruce: valor con distrbución uniforme $U(0.5, 1)$.
- Probabilidad de mutación: valor con distrbución uniforme $U(0.01, 0.2)$.
library(dplyr) set.seed(1234) n <- 100 tam_poblacion <- sample(10:200, size = n, replace = TRUE) prob_cruce <- runif(n = n, min = 0.5, max = 1) tam_torneo <- sample(2:10, size = n, replace =TRUE) prob_mutacion <- runif(n = n, min = 0.001, max = 0.2) grid <- data_frame(tam_poblacion, prob_cruce, tam_torneo, prob_mutacion) head(grid)
Cada combinación de parámetros se ha replicado $10$ veces para poder obtener datos sobre los que realizar análisis estadísticos.
La metodología que se ha llevado ha sido la de lanzar
Conclusiones
La resolución del Sudoku mediante el algoritmo genético implementado no parece ser la técnica más adecuada. En muchas de la mayoría de las combinaciones, el problema del Sudoku no estaba resuelto en las $5000$ primeras iteraciones, por lo que no se ha observado una eficiencia.
Una parte que se podría tratar de estudiar es la de si la representación de los individuos mediante un gen en cada celda con un número comprendido entre 1 y 9 es la más adecuada. Se propone como trabajo futuro el estudio de otras respresentaciones, como la representación con permutaciones. Es decir, cada individuo estaría representado por un genotipo $$ \begin{aligned} s = (&1, 2, \ldots, 9,\ &1, 2, \ldots, 9,\ &1, 2, \ldots, 9,\ &\ldots,\ &1, 2, \ldots, 9) \end{aligned} $$
Y cualquier individuo debería ser una permutación, lo que supondría elegir otro tipo de mutación que la random resetting elegida para esta implementación.
Referencias
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.