Estrellas Mágicas

miércoles, 5 de noviembre de 2008
En PC++ hay personajes curiosos. Uno de ellos es DjGera, que hace unos días propuso este problema:

Ubicar los números del 1 al 16 de forma que todas las columnas, filas y diagonales sumen 34.

Como buen ser humano mi primer instinto es mirar el problema fijo a ver si se me ocurre una solución. Pasados 10 minutos me empieza a dar hambre y recuerdo que para algo estudio Lic. en Computación.

Rápidamente, hay 16! = 2.1e13 formas de poner los números en el dibujo. Si obviamos las rotaciones y reflexiones de una misma solución, quedan apenas 653 837 184 000 formas posibles (seiscientos mil millones aproximadamente). Si bien las computadoras de hoy son rápidas, los problemas con números de esta magnitud todavía no son razonables para tratarlos "a lo bruto".

Solución con Programación Lineal

Al mirar el problema uno piensa bastante rápido (o no) en esas clases de álgebra lineal y empieza a imaginar matrices ralas y triangulaciones. Las restricciones del tipo "esta variable y esta otra deben sumar tanto" son efectivamente lineales y por lo tanto pueden representarse de forma matricial. Sin embargo, el problema como está planteado no es un sistema lineal convencional, por dos razones:
  • Las variables son enteras, no reales
  • Las variables deben ser todas distintas

Con estos añadidos, el problema se sale de la programación lineal y entra en la programación entera, un área que es computacionalmente mucho más difícil. Estos temas son objeto de la materia Investigación Operativa que es una de las gemas del Departamento de Computación de la UBA. Desgraciadamente, todavía no tuve oportunidad de cursar esa materia. Pensándolo un poco llegué a la conclusión de que probablemente esta no fuera la ruta más sabia, pero Martín, un poco más piola y sobre todo mucho más optimista, me mostró que estaba equivocado.

La técnica para resolver el problema usando programación entera es esencialmente:
  1. Plantear la matriz asociada a las restricciones lineales
  2. Triangularla con el método de eliminación gaussiana
  3. Tras este paso, queda una matriz singular de rango 8
  4. Fijar 8 de las 16 variables en valores enteros de 1 a 16
  5. Resolver las 8 variables restantes; si sus valores son enteros, están entre 1 y 16 y son diferentes a los 8 valores elegidos en el punto anterior, encontramos una solución.
  6. Volver al paso 4 para cada valuación distinta de las 8 variables a fijar

El algoritmo es entonces una combinación de programación lineal con fuerza bruta, pero programado en Mathematica tarda unas 10 horas en hallar todas las soluciones posibles, lo cual es muy respetable.

Solución con Algoritmo Genético

Con el rabo entre las piernas frente a la completitud y elegancia de la solución propuesta por Martín, implementé un algoritmo genético para resolver el problema. Sin demasiado trabajo el algoritmo encuentra una solución muy rápido (en unos 5 segundos con Psyco). Este tipo de algoritmos no deja de sorprenderme, puesto que si se observa el código las heurísticas de mutación y apareamiento son totalmente primitivas, y sin embargo el rendimiento del mismo es excelente.

#!/usr/bin/python
# estrella.py
# Resuelve la estrella magica de 16 puntas con un algoritmo genetico.

import random
import sys

tam_poblacion = 400
natalidad = 0.15
mutacion = 0.5
max_generacion = 500

# restricciones lineales
m = [[0,1,3,4],
     [4,5,7,8],
     [8,9,11,12],
     [12,13,15,0],
     [14,15,1,2],
     [2,3,5,6],
     [6,7,9,10],
     [10,11,13,14]]

# funcion de aptitud
def fitness(cand):
    fit = 0
    for eq in m:
        fit_local = -34
        for j in eq:
            fit_local += cand[j]
        fit += abs(fit_local)
    
    # esto fuerza que el primer numero sea el 16
    # y filtra las rotaciones de la misma solucion
    if cand[0] != 16:
        fit += 1000

    return fit

# funcion de apareamiento
def cruzar(c1, c2):
    l = [c1,c2]
    random.shuffle(l)
    c1, c2 = l
    c = c1[:8]
    for each in c2:
        if not each in c:
            c.append(each)
    return c

# funcion de mutacion
def mutar(cand):
    # "ejercito" al generador pseudoaleatorio
    # esto mejora bastante la velocidad de convergencia!
    i1 = random.randint(0,15)
    i2 = random.randint(0,15)
    
    i1 = random.randint(0,15)
    i2 = random.randint(0,15)
    
    i1 = random.randint(0,15)
    i2 = random.randint(0,15)

    c = cand[:]
    c[i1], c[i2] = c[i2], c[i1]
    
    return c

def avanzar_gen(pob):
    nacen = int(tam_poblacion * natalidad)
    mutan = int(tam_poblacion * mutacion)

    # unos se aparean
    for i in range(nacen):
        papa, mama = random.sample(poblacion, 2)
        poblacion.append(cruzar(papa, mama))

    # nacen unos mutados
    for i in range(mutan):
        deforme = random.sample(poblacion, 1)[0]
        poblacion.append(mutar(deforme))

    # ordeno por fitness
    poblacion.sort(lambda x,y: cmp(fitness(x), fitness(y)))
    
    # mato a los menos aptos
    while(len(poblacion) > tam_poblacion):
        poblacion.pop()
    

numeros = range(1,17)

# inicializador de poblacion
def init_pob():
    poblacion = []
    for i in range(tam_poblacion):
        cand = numeros[:]
        random.shuffle(cand)
        poblacion.append(cand)
    return poblacion


if __name__ == '__main__':
    while True:
        print
        print "Comenzando simulacion..."
        gen = 0
        poblacion = init_pob()
        min_fit = 1000
        while min_fit !=0:
            avanzar_gen(poblacion)
            gen += 1
            f = fitness(poblacion[0])
            
            if f < min_fit:
                print "Fitness: %s (gen %s)" % (f, gen)
            min_fit = f
            
            if min_fit == 0:
                print "Eureka!"
                print poblacion[0]
                sys.exit(0)
            if poblacion[0] == poblacion[len(poblacion)-1] or gen >= max_generacion:
                break
La ejecución se ve aproximadamente así (con un poco de suerte!):

gomo@ciabatta:~$ time python estrella.py

Comenzando simulacion...
Fitness: 30 (gen 1)
Fitness: 18 (gen 4)
Fitness: 16 (gen 10)
Fitness: 10 (gen 15)
Fitness: 8 (gen 30)
Fitness: 6 (gen 44)
Fitness: 4 (gen 47)
Fitness: 2 (gen 99)
Fitness: 0 (gen 237)
Eureka!
[16, 3, 11, 2, 13, 7, 14, 9, 5, 10, 1, 15, 4, 6, 12, 8]

real 0m2.449s
user 0m2.420s
sys 0m0.028s
Por su naturaleza, el tiempo que insume el algoritmo genético es muy variable y no encuentra siempre la misma solución. Por esta razón el programa hace varios intentos en lugar de perseverar con la misma simulación si observa que está tardando mucho en encontrar una solución.

Agregado el 10/12/2008
Solución con Backtracking

Tras discutir en un foro, considero que es mejor aclarar lo siguiente: los algoritmos propuestos arriba no son ni de cerca los más eficientes que se conocen para este problema. Un simple algoritmo de backtracking encuentra una solución de forma certera y bastante más veloz que lo que propongo arriba.

El interés del post era mostrar dos técnicas interesantes para resolver el problema, sin ánimos de indicar que eran particularmente aptas o buenas para la tarea. A modo de ilustración, adjunto el código del algoritmo de backtracking.

def buscarSolucion(actual, restantes):
    if restantes == []:
        return actual
    else:
        for r in restantes:
            nuevo_act = actual + [r]
            nuevo_rest = restantes[:]
            nuevo_rest.remove(r)
            
            s = None
            if esValidoParcial(nuevo_act):
                s = buscarSolucion(nuevo_act, nuevo_rest)
                
            if not s is None:
                return s
            

def esValidoParcial(perm):
    l = [[0,1,3,4],
         [4,5,7,8],
         [8,9,11,12],
         [12,13,15,0],
         [14,15,1,2],
         [2,3,5,6],
         [6,7,9,10],
         [10,11,13,14]]

   
    for rest in l:
        s = 0
        c = 0
        for k in rest:
            if len(perm) > k:
                s += perm[k]
                c += 1
        
        if c == 4 and s != 34:
            return False
        elif c < 4 and s >= 34:
            return False
    
    return True

if __name__ == "__main__":
    numeros = range(1,17)
    print buscarSolucion([], numeros)

2 comentarios:

Anónimo dijo...

Excelente explicacion!
Panchex

Anónimo dijo...
Este comentario ha sido eliminado por un administrador del blog.

Publicar un comentario