Tudo na linha....

Quando efectuamos experiências estatísticas podemos ter situações em que é o mesmo sujeito que é testado em diferentes situações ou não. Para cada uma das situações existem testes apropriados. No caso do SPSS cada linha de uma tabela corresponde aos diferentes resultados de um mesmo sujeito às diferentes condições. Isso implica que se tivermos sujeitos diferentes a testar as diferentes condições temos que colocar os resultados numa (longa ) coluna, acrescentado ainda uma segunda coluna de agrupamento. Nas experiências que fazemos, o que determina se estamos a fazer testes emparelhados ou não, é a natureza da população inicial em cada execução. Por exemplo, se estivermos a comparar o efeito de usar dois tipos diferentes de operadores de cruzamento e efectuarmos 30 execuções de cada condição, estaremos a efectuar testes emparelhados se em cada uma das 30 execuções as duas alternativas forem iniciadas com a mesma população. Acontece que, tipicamente, guardamos os resultados das diferentes execuções como se es testes fossem emparelhados. No caso de não serem temos a extenuante missão de passar tudo para uma coluna e introduzir a variável de agrupamento. Extenuante, a menos que tenhamos um programa que automatiza este processo.

def transform_file_spss(file_in, file_out):
    """From a file with matched values to a one without."""
    f_in = open(file_in)
    header_in = f_in.readline()
    body_in = f_in.readlines()
    f_in.close()
    
    split_header_in = header_in.split()
    number_col = len(split_header_in)
    number_lin = len(body_in)
    
    f_out = open(file_out,'w')
    header_out = 'Group\tValue\n'
    f_out.write(header_out)
    aux = [line[:-1].split() for line in body_in]
    # rearranje data
    body_out = []
    for i in range(number_col):
        # obtain column
        col = []
        for j in range(number_lin):
            col.append(aux[j][i])
        body_out.append(col)
    # write data
    for i in range(number_col):
        for j in range(number_lin):
            value = body_out[i][j]
            line_out = '%s\t%s\n' % (i,value)
            f_out.write(line_out)
    f_out.close() 

Nem mais!

Separações e outras complicações

Todos sabemos que os ficheiros existem sob diversos formatos. Os ficheiros de texto são muito simples e permitem guardar informação importante relativa à execução dos nossos algoritmos evolucionários. Os dados recolhidos são posteriormente tratados estatisticamente, sendo para isso necessário que o programa usado consiga importar correctamente o nosso ficheiro. Para tal temos que ter em atenção, entre outras coisas, dois aspectos: qual o símbolo usado para separar a parte inteira da parte decimal dos números e, qual o caracter separador dos diferentes valores. No que se refere ao SPSS (ou PASW) existe um comando que nos permite verificar qual o símbolo usado para separar as partes inteira e decimal (SHOW DECIMAL) e um modo de definir esse símbolo (SET DECIMAL {DOT, COMMA}). Basta abrir uma nova janela de sintaxe (File -> New -> Syntax), escrever e executar os comandos. Para facilitar a resolução da segunda questão podemos escrever um pequeno programa em Python que trata da questão.

def convert_file(file_in, sep_in,file_out,sep_out):
    """
    Change the separator of elements in a file.
    Assume ther existence of a header. Warning: the separator must be just that: a separator!
    """
    f_in = open(file_in)
    header_in = f_in.readline()
    body_in = f_in.readlines()
    f_in.close()
    
    f_out = open(file_out,'w')
    f_out.write(header_in)
    body_out = [line_in.replace(sep_in,sep_out) for line_in in body_in]
    f_out.writelines(body_out)
    f_out.close()

Simples, não é?

Danças Latinas

Um método de gerar uma população inicial é conhecido por latin hypercube. A ideia é simples: dividir o espaço em hipercubos todos iguais e gerar um indivíduo para cada hipercubo. O algoritmo apresentado nas aulas tinha um problema, pois fixava uma das dimensões, fazendo com que os valores escolhidos ao longo de uma das dimensões fosse sempre o mesmo. A figura ilustra o problema para o caso de duas dimensões.



Vamos por isso corrigir o código apresentado.

def latin_hypercube(domain, size):
    """One individual by hypercube."""
    numb_segments = size ** (1.0/len(domain))
    vectors = []
    for dom_i in domain:
        dist = (dom_i[1] - dom_i[0])/ numb_segments
        vectors.append(list(frange(dom_i[0],dom_i[1], dist))) 
    pairs = []
    for vector in vectors:
        pairs_from_vector = [[vector[i], vector[i+1]] for i in range(len(vector)-1)]
        pairs.append(pairs_from_vector)
    hyper_cubes = gen_cross_all(pairs)
    data_points = [generate_one(hc) for hc in hyper_cubes]
    return data_points  

Agora, depois de definirmos o número de segmentos, calculamos as coordenadas de todos os hipercubos, que armazenamos em hyper_cubes. De seguida fabricamos os pontos dentro de cada hipercubo, que armazenamos em data_points.

Cada ponto é gerado a partir de uma distribuição uniforme.

def generate_one(domain):
    return [random.uniform(dom_i[0],dom_i[1]) for dom_i in domain]

A geração das coordenadas para os hipercubos tem a sua dificuldade se queremos ter um programa genérico que funcione para uma dimensão qualquer e não apenas para duas dimensões.

A solução mais simples é recursiva.

def gen_cross_all(values):
    partial = [[elem] for elem in values[0]]
    remaining = values[1:]
    return gen_cross(partial, remaining)

def gen_cross(partial,remaining):
    """partial = list of points with dimension k.
    remaining = list of n-k dimensions."""
    if remaining == []:
        return partial
    else:
        partial = expand(partial,remaining[0])
        return gen_cross(partial,remaining[1:])
    

def expand(partial, dimension):
    new_partial = []
    for segment in dimension:
        new_partial.extend(cross(partial,segment))
    return new_partial

def cross(partial,segment):
    return [elem + [segment] for elem in partial]

A função gen_cross recebe uma lista parcial, com as coordenadas para um subconjunto k de dimensões, e a lista das dimensões restantes. Esta função é chamada por gen_cross_all, que converte a primeira dimensão numa primeira lista parcial e invoca depois a função que vai fazer todo o trabalho, isto é gen_cross.

Ao executar o código corrigido já conseguimos uma inicialização da população com as propriedades pretendidas.

Inimigos do infinito

Vimos na aula métodos de inicializar a população que permitem gerar indivíduos bem espalhado no espaço de procura. Um desse métodos baseia-se numa ideia muito simples: gerar um primeiro indivíduo aleatoriamente e, de seguida, ir acrescentando novos indivíduos de modo que cada um que é acrescentado mantém uma distância mínima a todos os outros. Relembremos o código.

def ssi(domain,size, min_dist):
    """not controlling max trials."""
    pop = [generate_one(domain)]
    for i in range(size-1):
        indiv = generate_one(domain)
        while not accept(indiv,pop,min_dist):
            indiv = generate_one(domain)
        pop.append(indiv)
    return pop

def accept(indiv,population, distance):
    for elem in population:
        if euclidean_distance(indiv,elem) < distance:
            return False
    return True

Como foi observado na aula, existe o risco de não conseguirmos gerar um novo indivíduo que satisfaça a condição de afastamento. A ser assim, o nosso programa nunca pára pois não conseguew abandonar o ciclo while. Vejamos como podemos remediar essa situação.

def ssi_b(domain_size, min_dist, max_reject):
    """Controlling max trials."""
    pop = [generate_one(domain)]
    for i in range(size-1):
        reject = 0
        indiv = generate_one(domain)
        while not accept(indiv,pop,min_dist):
            if reject == max_reject:
                return []
            else:
                reject += 1
                indiv = generate_one(domain)
        pop.append(indiv)
    return pop 

A ideia é trivial: fixamos um número máximo de tentativas. Se esse valor for ultrapassado o programa termina devolvendo uma população vazia. Claro que isto não resolve a questão essencial de precisarmos mesmo de uma população inicial! Por isso, teremos que lançar este programa o número de vezes necessário até que seja gerada uma população satisfazendo o pretendido. Cada vez que repetimos o programa, podemos, por exemplo, diminuir a distância que separa os indivíduos entre si, para aumentar a possibilidade de sucesso.

De alhos a bugalhos

Vamos falar de um problema... ou talvez não.

Existem ficheiros em formato tsp que como o nome indicam têm informação relevante para resolver o problema do caixeiro viajante. Eis um exemplo para o caso das cidades do Sahara Ocidental.

NAME : wi29
COMMENT : 29 locations in Western Sahara
COMMENT : Derived from National Imagery and Mapping Agency data
TYPE : TSP
DIMENSION : 29
EDGE_WEIGHT_TYPE : EUC_2D
NODE_COORD_SECTION
1 20833.3333 17100.0000
2 20900.0000 17066.6667
3 21300.0000 13016.6667
4 21600.0000 14150.0000
5 21600.0000 14966.6667
6 21600.0000 16500.0000
7 22183.3333 13133.3333
8 22583.3333 14300.0000
9 22683.3333 12716.6667
10 23616.6667 15866.6667
11 23700.0000 15933.3333
12 23883.3333 14533.3333
13 24166.6667 13250.0000
14 25149.1667 12365.8333
15 26133.3333 14500.0000
16 26150.0000 10550.0000
17 26283.3333 12766.6667
18 26433.3333 13433.3333
19 26550.0000 13850.0000
20 26733.3333 11683.3333
21 27026.1111 13051.9444
22 27096.1111 13415.8333
23 27153.6111 13203.3333
24 27166.6667 9833.3333
25 27233.3333 10450.0000
26 27233.3333 11783.3333
27 27266.6667 10383.3333
28 27433.3333 12400.0000
29 27462.5000 12992.2222
EOF

Depois de um cabeçalho aparecem as coordenadas das cidades precedidas de um identificador, no caso um inteiro com o número de ordem da cidade. Os algoritmos evolucionários que desenvolvemos para este problema assumem que o que conhecemos é a matriz das distâncias entre pares de cidades. Coloca-se por isso a necessidade de fazer um interface entre os dados guardados num ficheiro e as distâncias guardadas numa lista de listas. É o que é feito no código que se segue.

def le_coordenadas_tsp(ficheiro):
    """ Devolve a matriz das coordenadas a partir de ficheiro em formato tsp."""
    fich_in = open(ficheiro)
    tudo = fich_in.readlines()
    coordenadas = []
    for linha in tudo[7:-1]:
 n,x,y = linha[:-1].split()
 coordenadas.append((float(x),float(y)))
    fich_in.close()
    return coordenadas

def dicio_cidades(coordenadas):
    """ Cria um dicionário com as cidades e suas coordenadas."""
    dicio = {}
    for i, (x,y) in enumerate(coordenadas):
 dicio[i] = (x,y)
    return dicio


def fenotipo(genotipo,dic_cidades):
    """ Devolve o fenótipo associado."""
    fen = [dic_cidades[cidade] for cidade in genotipo]
    return fen

def evaluate(caminho):
    num_cidades = len(caminho)
    comp = 0
    for i in range(num_cidades):
 j = (i+1) % num_cidades
 comp  += distancia(caminho[i],caminho[j])
    return comp

def distancia(cid_i, cid_j):
    """ Distância Euclidiana."""
    x_i, y_i = cid_i
    x_j, y_j = cid_j
    dx = x_i - x_j
    dy = y_i - y_j
    dist = sqrt(dx**2 + dy**2)
    return dist

Na primeira parte a informação é extraída do ficheiro (le_coordenadas) e colocada num dicionário (dicio_cidades). De seguida podemos transformar uma permutação numa lista de coordenadas (fenotipo). A partir desse conhecimento podemos calcular as distância euclidiana entre qualquer par de cidades (evaluate e distancia). Nesta solução este código deve estar embebido no programa. Para ser tudo transformado de uma só vez, temos que fazer diferente:

def data_tsp_matrix(file_data):
    coordenadas = le_coordenadas_tsp(file_data)
    dic_cidades = dicio_cidades(coordenadas)
    city_matrix = [[0]* len(dic_cidades) for i in range(len(dic_cidades))]
    for i in range(len(dic_cidades)):
 for j in range(i+1,len(dic_cidades)):
     city_matrix[i][j] = city_matrix[j][i] = distancia(dic_cidades[i],dic_cidades[j])
    return city_matrix  

Este programa funciona para qualquer ficheiro em formato tsp. Admite a versão simétrica, isto é, ir de A a B custa o mesmo que ir de B a A.

Permuta-me as viagens

Em post anterior já falámos no problema do caixeiro viajante. Agora mostramos uma solução alternativa, recorrendo a permutações de inteiros. Os operadores de variação são a mutação por troca e o cruzamento por ordem. Sem mais explicações deixamos o código.

""" Basic Evolutionary Algorithms for the TSP problem.
Representation = integers permutation. """

__author__ = 'Ernesto Costa'
__date__ ='April 2012'

import math
import random
import copy
import operator

# ---
import matplotlib.pyplot as plt

# Main function

def ea_tsp_permutation(numb_cities,fitness,population_size,sel_parents,sel_survivals,mutation,prob_mutation,crossover, prob_crossover,num_generations):
    """
    numb_cities = number of cities
    fitness = to calculate the fitness : city matrix
    """
    population = random_population(numb_cities,population_size)
    population = evaluate_population(population,fitness)
    list_best_generation = [best_population(population)[1]]
    for i in range(num_generations):
        #print 'gener %d\n' % i
        aux_pop = copy.deepcopy(population)
        mates = sel_parents(aux_pop)
        offspring = crossover(mates, prob_crossover)
        offspring = mutation(offspring,prob_mutation)
        offspring = evaluate_offspring(offspring,fitness)
        population = sel_survivals(population,offspring)
        list_best_generation.append(best_population(population)[1])
    print 'FINAL BEST: ',list_best_generation[-1]
    return list_best_generation

# Generate population
def random_population(dim, pop_size):
    pop = []
    for i in range(pop_size):
 indiv = range(dim)
 random.shuffle(indiv)
 pop.append([indiv,0])
    return pop

def best_population(population):
    pop = copy.deepcopy(population)
    pop.sort(key= operator.itemgetter(1))
    return pop[0]
    
    
# Selection parents
def sel_parents(population):
    """Tournament of size 5."""
    new_pop = []
    for i in range(len(population)):
 indiv = tournament(population, 5)
 new_pop.append(indiv)
    return new_pop

def tournament(population, size):
    group = random.sample(population,size)
    group.sort(key=operator.itemgetter(1))
    return group[0]

# Mutation
def mutation(offspring, prob_mutation):
    new_off = [mutation_one(child,prob_mutation) for child in offspring]
    return new_off

def mutation_one(off,prob_mutation):
    """ Swap mutation."""
    muta = random.random()
    if muta < prob_mutation:
 cromo = off[0]
 i,j = random.sample(range(len(cromo)),2)
 cromo[i],cromo[j] = cromo[j],cromo[i]
 return [cromo,0]
    else:
 return off

# Crossover

def crossover(mates,prob_crossover):
    random.shuffle(mates)
    off = []
    for i in range(0,len(mates)-1,2):
 off.extend(order_cross(mates[i][0],mates[i+1][0],prob_crossover))
    return off

# OX - order crossover
def order_cross(cr_1,cr_2,prob_cross):
 size = len(cr_1)
 value = random.random()
 if value < prob_cross:
     pc= random.sample(range(size),2)
     pc.sort()
     pc1,pc2 = pc
     f1 = [None] * size
     f2 = [None] * size
     f1[pc1:pc2+1] = cr_2[pc1:pc2+1]
     f2[pc1:pc2+1] = cr_1[pc1:pc2+1]
     indexes = range(pc2+1,size) + range(pc1)
     j = k = pc2 + 1
     for i in indexes:  
  while cr_1[j%size] in f1:
      j += 1
  f1[i%size] = cr_1[j%size]
  j += 1
  while cr_2[k%size] in f2:
      k += 1
  f2[i%size] = cr_2[k%size]
  k += 1
     return [[f1,0],[f2,0]]
 else:
     return [[cr_1,0],[cr_2,0]]

 
# Selection of the survivals
def sel_survivals(population,offspring):
    """ mu + lambda."""
    all_pop = population + offspring
    all_pop.sort(key=operator.itemgetter(1))
    return all_pop[:len(population)]
 
def sel_survivals_elite(population,offspring):
    """ Elitism."""
    population.sort(key=operator.itemgetter(1))
    offspring.sort(key=operator.itemgetter(1), reverse=True)
    elite = int(0.05 * len(population))
    all_pop = population[:elite] + offspring[elite:]
    random.shuffle(all_pop)
    return all_pop


# Evaluation of the fitness
def evaluate_population(population,fitness):
    population = [[individual[0],evaluate(individual[0], fitness)] for individual in population]
    return population

def evaluate_offspring(offspring,fitness):
    return evaluate_population(offspring,fitness)
 
def evaluate(tour, city_matrix):
    numb_cities = len(tour)
    cost = 0
    for i in range(numb_cities):
 j = (i+1) % numb_cities
 cost  += city_matrix[tour[i]][tour[j]]
    return cost


# Create data if needed
def create_distance_matrix(size):
    matrix = [ [0] * size for i in range(size)]
    for i in range(size):
 for j in range(i+1,size):
     value = random.randint(500,1500)
     matrix[i][j] = matrix[j][i] = value
    return matrix

# Simple visualization
   
def visualize_fitness(list_values):
    x = range(1,len(list_values)+1)
    y = list_values
    plt.title('Best Fitness')
    plt.xlabel('Generation')
    plt.ylabel('Fitness')
    plt.plot(x,y,'r-o')
    #plt.savefig('one_and_one.png')
    plt.show()
 
      
def main(num_runs, algorithm,problem):
    all_best = []
    for i in range(num_runs):
        best = algorithm(problem[0],problem[1],150, sel_parents, sel_survivals_elite,mutation,\
                          0.01, crossover, 0.8,5000)
        all_best.append(best)
    transformed = zip(*all_best)
    best_best = [max(gener) for gener in transformed]
    average_best = [sum(gener) / float(num_runs) for gener in transformed]
    visualize_fitness(average_best)  

Chaves aleatórias

O problema do caixeiro viajante (TSP) é nosso velho conhecido: alguém pretende visitar um certo número de cidades, uma e uma só vez, regressando à cidade de origem. Quer minimizar o custo da viagem, entendido aqui pela menor distância percorrida. Como em qualquer problema envolvendo algoritmos evolucionários a primeira questão que se coloca é a da representação das soluções candidatas. A solução mais evidfewnte é uma permutação de inteiros, que nos obriga a escolher posteriormente com cuidado os operadores de variação. Mas existem alternativas. Uma recorre às chamadas chaves aleatórias (random keys). Trata-se de vectores de reais, entre 0 e 1 que permitem representar permutações. A vantagem é a de que se podem usar os operadores de variação simples para reais. Claro que necessitamos de estabelecer uma relação entre um vector de reais e uma permutação. O código para o fazer é trivial:

def decode_rk(vector):
    aux = deepcopy(vector)
    aux.sort()
    permutation =[ vector.index(elem) + 1 for elem in aux]
    return permutation

Ordenamos o vector e procuramos a sua posição no vector original!

O algoritmo para o TSP segue o modelo habitual de um algoritmo evolucionário.

def ea_tsp_rk(numb_cities,fitness,population_size,sel_parents,sel_survivals,mutation,
prob_mutation,crossover, prob_crossover,num_generations):
    """
    numb_cities = number of cities
    fitness = to calculate the fitness : city matrix
    """
    population = random_population(numb_cities,population_size)
    population = evaluate_population(population,fitness)
    list_best_generation = [best_population(population)[1]]
    for i in range(num_generations):
        print 'gener %d\n' % i
        aux_pop = copy.deepcopy(population)
        mates = sel_parents(aux_pop)
        offspring = crossover(mates, prob_crossover)
        offspring = mutation(offspring,prob_mutation)
        offspring = evaluate_offspring(offspring,fitness)
        population = sel_survivals(population,offspring)
        list_best_generation.append(best_population(population)[1])
    print 'FINAL BEST: ',list_best_generation[-1]
    return list_best_generation

O algoritmo utilizado recorre a uma selecção por torneio, mutação gaussiana, cruzamento uniforme e selecção de sobreviventes de estado estável.

Selecção por torneio de tamanho 5

def sel_parents(population):
    """Tournament of size 5."""
    new_pop = []
    for i in range(len(population)):
 indiv = tournament(population, 5)
 new_pop.append(indiv)
    return new_pop

def sel_tournament(population,size):
    """Tournament of size."""
    new_pop = []
    for i in range(len(population)):
 indiv = tournament(population, size)
 new_pop.append(indiv)
    return new_pop

def tournament(population, size):
    """Minimization problem."""
    group = random.sample(population,size)
    group.sort(key=operator.itemgetter(1))
    return group[0]

Mutação Gaussiana

def mutation(offspring, prob_mutation):
    new_off = [mutation_one(child,prob_mutation) for child in offspring]
    return new_off

def mutation_one(off,prob_mutation):
    """ Gaussian mutation."""
    muta = random.random()
    if muta < prob_mutation:
 cromo = off[0]
 index = random.randint(0,len(cromo)-1)
 value = cromo[index]
 delta = random.gauss(0,1.0)
 while not in_bounds(value + delta):
     delta = random.gauss(0,1.0)
 cromo[index] = value + delta
 return [cromo,0]
    else:
 return off

def in_bounds(value):
    return (0.0 <= value <= 1.0)

Cruzamento Uniforme

def uniform_crossover(mates,prob_crossover):
    random.shuffle(mates)
    off = []
    for i in range(0,len(mates)-1,2):
 off.extend(uni_cross(mates[i][0],mates[i+1][0],prob_crossover))
    return off

def uni_cross(cromo_1, cromo_2,prob_cross):
 value = random.random()
 if value < prob_cross:
  f1=[]
  f2=[]
  for i in range(0,len(cromo_1)):
   if random.random() < 0.5:
    f1.append(cromo_1[i])
    f2.append(cromo_2[i])
   else:
    f1.append(cromo_2[i])
    f2.append(cromo_1[i])

  return [[f1,0],[f2,0]]
 else:
  return [[cromo_1,0],[cromo_2,0]]

Selecção Estado Estável

def sel_survivals(population,offspring):
    """ mu + lambda."""
    all_pop = population + offspring
    all_pop.sort(key=operator.itemgetter(1))
    return all_pop[:len(population)]

Avaliação

A avaliação é feita no fenótipo, isto é, numa permutação de inteiros, usando uma matriz de distâncias.

def evaluate_population(population,fitness):
    population = [[individual[0],evaluate(fenotype(individual[0]), fitness)] for individual in population]
    return population

def evaluate_offspring(offspring,fitness):
    return evaluate_population(offspring,fitness)

def fenotype(genotype):
    """From a random key to a permutation."""
    aux = copy.deepcopy(genotype)
    aux.sort()
    permutation =[ genotype.index(elem)  for elem in aux]
    return permutation    

def evaluate(tour, city_matrix):
    numb_cities = len(tour)
    cost = 0
    for i in range(numb_cities):
 j = (i+1) % numb_cities
 cost  += city_matrix[tour[i]][tour[j]]
    return cost

Testar

Pode testar o algoritmo gerando a sua própria matriz de distâncias:

def create_distance_matrix(size):
    matrix = [ [0] * size for i in range(size)]
    for i in range(size):
 for j in range(i+1,size):
     value = random.randint(500,1500)
     matrix[i][j] = matrix[j][i] = value
    return matrix

Visualizar

Para ver o resultado:

def visualize_fitness(list_values):
    x = range(1,len(list_values)+1)
    y = list_values
    plt.title('Best Fitness')
    plt.xlabel('Generation')
    plt.ylabel('Fitness')
    plt.plot(x,y,'r-o')
    #plt.savefig('one_and_one.png')
    plt.show()

Executar

Para executar (cuidasdo que pode demorar muito tempo....)

def main_all(num_runs, algorithm,problem):
    all_best = []
    for i in range(num_runs):
        best = algorithm(problem[0],problem[1],50, sel_parents, sel_survivals,mutation,\
                          0.1, uniform_crossover, 0.7,1000)
        all_best.append(best)
    transformed = zip(*all_best)
    best_best = [max(gener) for gener in transformed]
    average_best = [sum(gener) / float(num_runs) for gener in transformed]
    visualize_fitness(average_best) 

Exemplo

No nosso caso usámos a matriz das distâncias das 29 cidades do Sahara Ocidental (ver ficheiro wi29.tsp na InforEstudante ou procure pelo nome usando um motor de busca. este ficheiro guarda as coordenadas das cidades que têm que ser transformadas em distâncias). O resultado óptimo é 27603.

A solução óptima.


Uma execução.

Lei da Morte

Depois de construirmos uma nova geração, quem de entre pais e filhos, sobrevive? Alguns métodos.

Geracional

Nadas mais simples e fácil do que sobreviverem apenas os filhos. O único problema é perderem-se pais de melhor qualidade. Mas não se pode ter tudo...

def survivors_generational(parents,offspring):
    return offspring

Elitismo

Mas se quisermos garantir que não se perdem os melhores pais podemos garantir que eles passam sempre, sem alteração, para a geração seguinte. O número dos elementos dos pais que passam para a geração seguinte, o tamanho da élite, pode ser controlado. No limite basta apenas manter o melhor.

def survivors_elitism(parents,offspring,elite):
    """ Sorted. """
    size = len(parents)
    comp_elite = int(size* elite)
    new_population = parents[:comp_elite] + offspring[:size - comp_elite]
    return new_population

Estado Estável
Podemos colocar as duas gerações em competição pela sobrevivência, cabendo a decisão à qualidade de cada indivíduo.

def survivors_mu_plus_lambda(parents,offspring):
    """Minimizing."""
    size = len(parents)
    new_population = parents + offspring
    new_population.sort(key=itemgetter(1)) 
    return new_population[:size]

Não é difícil imaginar pequenas variantes a estes métodos.

Ser Selectivo

Escolher os melhores para reprodução com variação é um elemento chave dos algoritmos evolucionários. Tem as suas nuances, pois demasiada pressão selectiva pode causara convergência prematura da população para um máximo local. Vários mecanismos de selecção foram propostos ao longo do tempo. Vejamos alguns.

Torneio

Simples de implementar, procurando contrariar a pressão selectiva ao escolher de modo aleatório os elementos que vão participar no torneio.

def tournament_sel(population, numb,size):
    mate_pool=[]
    for i in range(numb):
        indiv = tournament(population,size)
        mate_pool.append(indiv)
    return indiv

def tournament(population, numb_offspring, t_size):
    return [one_tournament(population,t_size) for count in range(numb_offspring)]

def one_tournament(population,size):
    """Maximization. Problem.Deterministic"""
    pool = random.sample(population, size)
    pool.sort(key=itemgetter(1),reverse = True)
    return pool[0]

Trata-se de um torneio determinístico, significando isso que o melhor é sempre o escolhido. Existe uma versão estocástica que usa um processo semelhante ao da recristalização simulada, podendo ser escolhido, com uma dada probabilidade, um indivíduo de menor qualidade.


Ordem
Outro algoritmo que procura contrariar a pressão selectiva é o algoritmo de selecção por ordem. Ordenam-se os indivíduos em função da sua qualidade mas escolha depende da posição e não do valor do mérito.

def rank(population, numb,s):
    """linear version."""
    pop_size = len(population)
    pop = population[:]
    pop.sort(key=itemgetter(1))
    probabilities = [((2 - s)/pop_size) + ((2*i*(s -1))/float((pop_size * (pop_size - 1))))
                     for i in range(pop_size)]
    mate_pool = []
    for i in range(numb):
        value = random.uniform(0,1)
        index = 0
        total = probabilities[index]
        while total < value:
            index += 1
            total += probabilities[index]
        mate_pool.append(population[index])
    return mate_pool

Roleta

O algoritmo clássico. A escolha de um indíviduo é directamente proporcional ao seu mérito. Tem problemas de diminuição da diversidade da população ao longo das gerações e, por isso, de convergência prematura. Para seleccionar n indivíduos simula-se a rotação de uma roleta igual número de vezes.

def roulette_wheel(population, numb):
    """ Select numb individuals from the population
    according with their relative fitness. MAX. """
    pop = population[:]
    pop.sort(key=itemgetter(1))
    total_fitness = sum([indiv[1] for indiv in pop])
    mate_pool = []
    for i in range(numb):
        value = random.uniform(0,1)
        index = 0
        total = pop[index][1]/ float(total_fitness)
        while total < value:
            index += 1
            total += pop[index][1]/ float(total_fitness)
        mate_pool.append(pop[index])
    return mate_pool

Amostragem Universal Estocástica

Este método permite com uma rotação apenas da roleta escolher todos os elementos. É preferível ao método da roleta clássico sempre que queremos obter mais do que um elemento. Se tivermos que produzir n descendentes escolhe-se um ponto aleatório entre [0, 1/n]. Depois marcam-se posições igualmente espaçadas de 1/n. A roleta está dividida de acordo com o mérito dos indivíduos. Roda-se uma vez e escolhem-se os indivíduos em frente de vada uma das n marcas.

def sus(population,numb):
    """ Stochastic Universal Sampling."""
    pop = population[:]
    pop.sort(key=itemgetter(1))
    total_fitness = sum([indiv[1] for indiv in pop])
    value = random.uniform(0,1.0/numb)
    pointers = [ value + i * (1.0/numb) for i in range(numb)]
    mate_pool = []
    for j in range(numb):
        val = pointers[j]
        index = 0
        total =pop[index][1]/float(total_fitness)
        while total < val:
            index += 1
            total += pop[index][1]/float(total_fitness)
        mate_pool.append(pop[index])
    return mate_pool

Cruzamentos perigosos

O operador de cruzamento altera os indivíduos sem no entanto introduzir valores novos. Comecemos por por operadores genéricos, tipicamente usados com cadeias binárias e com vectores de reais. Se os pais forem válidos os filhos também o serão.

Um ponto de corte
Define um ponto aleatório e cruza.

def one_point_cross(cromo_1, cromo_2,prob_cross):
 value = random.random()
 if value < prob_cross:
  pos = random.randint(0,len(cromo_1))
  f1 = cromo_1[0:pos] + cromo_2[pos:]
  f2 = cromo_2[0:pos] + cromo_1[pos:]
  return [f1,f2]
 else:
  return [cromo_1,cromo_2]

Dois pontos de corte
Define dois pontos aleatórios e cruza.

def two_points_cross(cromo_1, cromo_2,prob_cross):
 print cromo_1, cromo_2
 value = random.random()
 if value < prob_cross:
  pc= random.sample(xrange(len(cromo_1)),2)
  pc.sort()
  pc1,pc2 = pc
  f1= cromo_1[:pc1] + cromo_2[pc1:pc2] + cromo_1[pc2:]
  f2= cromo_2[:pc1] + cromo_1[pc1:pc2] + cromo_2[pc2:]
  return [f1,f2]
 else:
  return [cromo_1,cromo_2]

Uniforme
Cada filho tem igual probabilidade de herdar de um dos pais.

def uniform_cross(cromo_1, cromo_2,prob_cross):
 value = random.random()
 if value < prob_cross:
  f1=[]
  f2=[]
  for i in range(0,len(cromo_1)):
   if random.random() < 0.5:
    f1.append(cromo_1[i])
    f2.append(cromo_2[i])
   else:
    f1.append(cromo_2[i])
    f2.append(cromo_1[i])

  return [f1,f2]
 else:
  return [cromo_1,cromo_2]

Qunado se trata de permutações temos que ter cuidados acrescidos para que do cruzamento resultem indivíduos inviáveis. Claro que, em muitas situações, a produção de indivíduos inviáveis é inevitável, e, em situações precisas, até pode ser desejável. Existem vários métodos para lidar com os indivíduos inviáveis que vão desde a sua penalização até à sua reparação. Mas regressemos aos operadores.

Ordem

Escolhem-se os mesmos dois pontos de corte nos dos pais P1 e P2. Trocam-se as partes do meio (as partes entre os dois pontos de corte: A parte do meio de P2 vai para F1 e a de P1 para F2. Em cada um dos filhos procura-se manter pela mesma ordem relativa, a partir do segundo ponto de corte os elementos do pai com o mesmo número de ordem. Por exemplo, No F1, que herda a parte do meio de P2, procura-se manter na mesma ordem os elementos de P1.

def order_cross(cr_1,cr_2,prob_cross):
 size = len(cr_1)
 value = random.random()
 if value < prob_cross:
  pc= random.sample(range(size),2)
  pc.sort()
  pc1,pc2 = pc
  f1 = [None] * size
  f2 = [None] * size
  f1[pc1:pc2+1] = cr_2[pc1:pc2+1]
  f2[pc1:pc2+1] = cr_1[pc1:pc2+1]
  indexes = range(pc2+1,size) + range(pc1)
  j = k = pc2 + 1
  for i in indexes:  
   while cr_1[j%size] in f1:
    j += 1
   f1[i%size] = cr_1[j%size]
   j += 1
   while cr_2[k%size] in f2:
    k += 1
   f2[i%size] = cr_2[k%size]
   k += 1
  return [f1,f2]
 else:
  return [cr_1,cr_2] 
  

Partial Mapped

Este operador começa também por determinar dois pontos de corte e trocar as partes do meio. De seguida procura manter o máximo de elementos fora da parte do meio nas suas posições iniciais. Nos casos em que tal não é possível, usa a correspondência desses elementos com os da outra parte do meio.





Não se assuste com o código.....

def pmx_cross(cromo_1,cromo_2,prob_cross):
 size = len(cromo_1)
 value = random.random()
 if value < prob_cross:
  pc= random.sample(xrange(size),2)
  pc.sort()
  pc1,pc2 = pc
  f1 = [None] * size
  f2 = [None] * size
  f1[pc1:pc2+1] = cromo_1[pc1:pc2+1]
  f2[pc1:pc2+1] = cromo_2[pc1:pc2+1]
  # primeiro filho
  # parte do meio
  for j in range(pc1,pc2+1):
   if cromo_2[j] not in f1:
    pos_2 = j
    g_j_2 = cromo_2[pos_2]
    g_f1 = f1[pos_2]
    index_2 = cromo_2.index(g_f1)
    while f1[index_2] != None:
     index_2 = cromo_2.index(f1[index_2])
    f1[index_2] = g_j_2
  # restantes
  for k in range(size):
   if f1[k] == None:
    f1[k] = cromo_2[k]
  # segundo filho 
  # parte do meio
  for j in range(pc1,pc2+1):
   if cromo_1[j] not in f2:
    pos_1 = j
    g_j_1 = cromo_1[pos_1]
    g_f2 = f2[pos_1]
    index_1 = cromo_1.index(g_f2)
    while f2[index_1] != None:
     index_1 = cromo_1.index(f2[index_1])
    f2[index_1] = g_j_1
  # parte restante
  for k in range(size):
   if f2[k] == None:
    f2[k] = cromo_1[k]    
  return [f1,f2]
 else:
  return [cromo_1,cromo_2]

Outros operadores

Existem muitos mais operadores de cruzamento. Uma vez mais, em muitos casos esses operadores estão ligados a problemas concretos. Deixamos aqui exemplos simples.

Aritmético

Cria dois descendentes baseado numa combinação linear dos seus pais.

def arithmetic(cromo_1,cromo_2,coef):
 off_1 = [ coef * cromo_1[index] + (1 - coef) * cromo_2[index] for index in range(len(cromo_1))]
 off_2 = [ coef * cromo_2[index] + (1 - coef) * cromo_1[index] for index in range(len(cromo_1))]
 return [off_1,off_2]

Heurístico

Cria apenas um descendente enviesado pelo valor do melhor dos pais.

def heuristic(cromo_1,cromo_2,coef,fit_func):
 """Just one offspring!."""
 fit_cr_1 = fitness(fit_func,cromo_1)
 fit_cr_2 = fitness(fit_func,cromo_2)
 if fit_cr_1 > fit_cr_2:
  off = [ coef* (cromo_1[index] - cromo_2[index]) + cromo_1[index] for index in range(len(cromo_1))]
 else:
  off = [ coef* (cromo_2[index] - cromo_1[index]) + cromo_2[index] for index in range(len(cromo_1))]
 return off

Mutemo-nos

Nota prévia: nos algoritmos apresentados as alterações são permanentes. Deve garantir que esta questão não coloca problemas ao código que os usa.


Comecemos pelo mais simples: mutação de cadeias binárias.

def muta_bin(cromo,prob_muta, muta_func):
    for i in range(len(indiv)):
        cromo[i] = muta_func(cromo[i],prob_muta)
    return cromo

# binário
def muta_bin_gene(gene, prob_muta):
    g = gene
    value = random.random()
    if value < prob_muta:
        g ^= 1
    return g

Não há muito a dizer: cada gene tem uma certa probabilidade de ser mutado. O valor da probabilidade de mutação é um parâmetro. Há quem use um valor baseado na regra de que que essa probabilidade deve ser o inverso do comprimento do indivíduo, garantindo assim que, em média, um gene é mutado.

Passemos ao caso de permutações de inteiros. Falámos nas aulas no operador que escolhe duas posições e troca os respectivos valores:

def muta_perm_swap(cromo, prob_muta):
    value = random.random()
    if value < prob_muta:
        index = random.sample(xrange(len(cromo)),2)
        index.sort()
        i1,i2 = index
        cromo[i1],cromo[i2] = cromo[i2], cromo[i1]
    return cromo

Mas existem alternativas, isto é, operadores que quando aplicados transformam uma permutação noutra permutação.

Escolhe dois pontos e baralha o conteúdo do vector entre os dois pontos:

def muta_perm_scramble(cromo,prob_muta):
    value = random.random()
    if value < prob_muta:
        index = random.sample(xrange(len(cromo)),2)
        index.sort()
        i1,i2 = index
        scramble = cromo[i1:i2+1]
        random.shuffle(scramble) 
        cromo = cromo[:i1] + scramble + cromo[i2+1:]
    return cromo

Escolhe dois pontos e inverte o conteúdo entre os dois pontos:

def muta_perm_inversion(cromo,prob_muta):
    value = random.random()
    if value < prob_muta:
        index = random.sample(xrange(len(cromo)),2)
        index.sort()
        i1,i2 = index
        inverte = []
        for elem in cromo[i1:i2+1]:
              inverte = [elem] + inverte
        cromo = cromo[:i1] + inverte + cromo[i2+1:]
    return cromo

Selecionar dois pontos. Insere o valor do gene na segunda posição no primeiro e desloca desde a primeira até à segunda uma posição para a direita:

def muta_perm_insertion(cromo, prob_muta):
    value = random.random()
    if value < prob_muta:
        index = random.sample(xrange(len(cromo)),2)
        index.sort()
        i1,i2 = index
        gene = cromo[i2]
        for i in range(i2,i1,-1):
          cromo[i] = cromo[i-1]
        cromo[i1+1] = gene
    return cromo

Finalmente, vectores de reais. Cado da mutação gaussiana.

def muta_reals(icromo, prob_muta, domain, sigma):
    for i in range(len(cromo)):
         cromo[i] = muta_reals_gene(cromo[i],prob_muta, domain[i], sigma[i])
    return cromo

Aqui verificamos se o novo valor excede o intervalo possível e, em caso afirmativo, aproxima esse valor do extremo mais próximo:

def muta_reals_gene_b(gene,prob_muta, domain_gene, sigma_gene):
    value = random.random()
    if value < prob_muta:
        muta_value = random.gauss(0,sigma_gene)
        new_gene = gene + muta_value
        if new_gene < domain_gene[0]:
            new_gene = domain_gene[0]
        elif new_gene > domain_gene[1]:
           new_gene = domain_gene[1]
    return new_gene

Enquanto o novo valor não for válido continua a gerar valores:

 
def muta_reals_gene(gene,prob_muta, domain_gene, sigma_gene):
    value = random.random()
    if value < prob_muta:
       muta_value = random.gauss(0,sigma_gene)
       new_gene = gene + muta_value
       while not (domain_gene[0] <= new_gene <= domain_gene[1]):
           muta_value = random.gauss(0,sigma_gene)
       new_gene = gene + muta_value             
    return new_gene

Aspectos de Projecto

Projectar um algoritmo evolucionário para um projecto concreto envolve aspectos de arte e de ciência (ou de design e de engenharia). Existem duas questões centrais a serem tratadas no início: que representação, que operadores de variação. Depois outros aspectos são relevantes: população inicial, métodos de selecção dos pais e dos sobreviventes, critério de paragem, os valores dos parâmetros tamanho da população, probabilidade de mutação e de cruzamento. É sobre alguns destes aspectos que os próximos posts vão discorrer. Para já estaremos limitados a representações binárias, permutações de inteiros, vectores de números reais.

Visões

Numa aula anterior andámos com o problema de desenvolver um algoritmo simples para reconhecer padrões. No caso analisado os padrões eram dígitos. O padrão era codificado como uma lista de listas. Por exemplo, o número um era codificado como:

[[0,0,1,0,0,],[0,1,1,0,0],[0,0,1,0,0],[0,0,1,0,0],[0,0,1,0,0]]

Vamos exemplificar como podemos mostrar o número como se de uma imagem se tratasse. para isso vamos pedir auxílio ao módulo pygame (pygame.org).

E primeiro lugar definimos uma classe Grid que representa os nossos padrões.

class Grid(object):
    def __init__(self, width,height):
        self.width = width
        self.height = height
        self.grid = [[0]* width for i in range(height)]
        
    def initialize(self):
        for x in range(self.width):
            for y in range(self.height):
                self.grid[y][x] = random.randint(0,1)
                
    def set_cell(self, x,y, value):
        self.grid[y][x] = value
                
    def set_grid(self,figure):
        for x in range(self.width):
            for y in range(self.height):
                self.set_cell(x,y,figure[y][x])
                
    def show_grid(self):
        for y in range(self.height):
            print self.grid[y]

Feito isto usamos o pygame para desenhar o padrão. A lógica é a de fazer corresponder um elemento da lista de listas a um quadrado na imagem com um número de pixeis configurável (ver imagem).





Os guiões (scripts) em pygame seguem a metodologia IDEA (Initialize, Display, Entities, Action). Passemos ao código.

def draw(figure):
    # Initialize
    WHITE = (255,255,255)
    BLACK = (0,0,0)
    SIZE_X = len(figure)
    SIZE_Y = len(figure[0])
    DISP_SIZE = 50
    # Display
    pygame.init()
    screen =  pygame.display.set_mode((SIZE_X * DISP_SIZE, SIZE_Y * DISP_SIZE))
    screen.fill(BLACK)
    # Entities
    figure_grid = Grid(SIZE_X, SIZE_Y)
    figure_grid.set_grid(figure)    
    # Action
    for x in range(figure_grid.width):
        for y in range(figure_grid.height):
            if figure_grid.grid[y][x] == 1:
                pygame.draw.rect(screen,WHITE, pygame.Rect(x*DISP_SIZE,y*DISP_SIZE, DISP_SIZE,DISP_SIZE))
                
    pygame.display.flip()

Executando o código:

if __name__ == '__main__':
    figure = [[0,0,0,1,1,0,0,0],[0,0,1,0,0,1,0,0],[0,0,0,0,0,1,0,0],[0,0,0,0,0,1,0,0],[0,0,0,0,1,0,0,0],[0,0,0,1,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,1,1,1,1,0,0]]
    draw(figure)

Obtemos o lindo padrão:

Uma imagem vale mais do que mil palavras

Em computação de inspiração biológicavisualizar é bastante importante. No passado já dissemos como podemos fazer o gráfico de funções unidimensionais, recorrendo ao módulo matplotlib. Foi o caso, por exemplo, da evolução da qualidade do melhor indivíduo da população ao longo das gerações. Mas também podemos ter interesse em ver a função objectivo associada a um problema. Vamos ver como podemos continuar a usar o módulo matlotlib para o fazer. Vejamos como exemplo a função Sphere para duas variáveis independentes:

import random 
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np


def show_3d(domain,function,title):
    # preparing
    fig = plt.figure()
    fig.suptitle(title, fontsize=18, fontweight='bold')
    ax = Axes3D(fig)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    # data
    X = np.arange(domain[0], domain[1], 0.15)
    Y = np.arange(domain[0], domain[1], 0.15)
    X, Y = np.meshgrid(X, Y)
    Z = function(X,Y)
    # rendering
    ax.plot_surface(X, Y, Z, rstride=1,cstride=1,cmap=cm.jet)
    ax.legend()
    plt.show()

# -- Problem
DOMAIN_DJ_1 = [-5,5]
    
def dejong_1_a(x,y):
    return x ** 2 + y ** 2
if __name__ == '__main__':
    show_3d(DOMAIN_DJ_1,dejong_1_a, 'Sphere')

Ao executar o código obtemos a figura retratada na imagem. Podemos verificar a existência de 4 máximos nas extremidades do domínio.




Em termos de código, e em alternativa, podemos recorrer a linspace em vez de arange. A diferença entre ambas, é que no caso da primeira dizemos quantos pontos queremos igualmente espaçados entre os limites inferior e superior.

def show_3d_b(domain,function,title):
    # preparing
    fig = plt.figure()
    fig.suptitle(title, fontsize=18, fontweight='bold')
    ax = Axes3D(fig)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')    
    # data
    X = np.linspace(-5,5,50)
    Y = np.linspace(-5,5,50)
    X, Y = np.meshgrid(X, Y)
    Z = function(X,Y)
    # rendering
    ax.plot_surface(X, Y, Z, rstride=1,cstride=1,cmap=cm.jet)
    ax.legend()
    plt.show()

Na aula propus duas variantes à função Sphere. Uma primeira em que alteramos a sua segunda componente, ficando agora a função a ser:

def dejong_1_b(x,y):
    return x ** 2 + y ** 3

Visualmente fica:


Uma alteração aparentemente menor modifica bastante o resultado final.
Uma segunda em que introduzimos ruído gaussiano:

def dejong_1_c(x,y):
    noise = np.array([random.gauss(0,3) for i in range(len(x))])
    return (x ** 2 + y ** 2) + noise

De novo mostramos a imagem. Neste caso o valor do ruído acrescentado altera a imagem obtida.

ART DECO

Os algoritmos evolucionários não se destinam a resolver todo o tipo de problemas. Estão vocacionados para resolver problemas que ou não têm solução analítica ou que, tendo, o espaço de soluções candidatas tem uma dimensão muito elevada, que torna intratatável computacionalmente o uso de técnicas exactas. Por isso, perante um problema, a primeira questão a resolver é a de se saber se os AE são ou não adequados. Se a resposta for afirmativa, a questão seguinte que se coloca é a de saber que instância concreta iremos usar. Para responder, recordemos a versão genérica de um AE.



Resulta claro que existem várias componentes do algoritmos que precisamos definir:

1- Como definir a população inicial
2- Como seleccionar os pais para reprodução
3- Como seleccionar de entre pais e filhos os sobreviventes para a geração seguinte
4- Que operadores genéticos vamos usar para alterar os pais que foram seleccionados
5- Como avaliar a qualidade das soluções candidatas
6- Como determinar o fim do processo evolutivo
7- Que parâmetros usar (e.g., tamanho da população, probabilidades dos operadores de variação)

Mas antes de mais temos que escolher uma representação para os indivíduos (cromossomas) da nossa população. Esta é, aliás, a questão primeira a resolver e, claramente, a mais importante de todas. As representações são essencialmente de dois tipos: lineares e não lineares. Podem ainda ser detalhadas: vectores binários, de inteiros, de reais, árvores, grafos, etc.. Uma boa representação deve ter quatro características:


1- económica: reduzir a dimensão do espaço de procura;
2- eficaz e deficiente: permitir encontrar a solução e encontrá-la em tempo curto (em associação com os operadores de variação);
3- completa: deve permitir representar todas as soluções;
4- ligada: ligar qualquer par de soluções (em associação com os operadores de variação)

A partir do momento que representação e operadores estão escolhidos, temos que ter em atenção as consequências, em particular saber distinguir:

- O espaço de procura conceptual;
- O espaço de procura real, consequência da representação;
- O espaço de procura efectivo, consequência dos operadores

Tipicamente estas três situações encontram-se numa situação de inclusão ou identidade. A figura abaixo tenta ilustrar a situação.





Exemplos

Vejamos exemplos concretos.

(1) Comecemos pelo problema OneMax. Neste exemplo o conjunto de todas as cadeias binárias é o espaço conceptual. Usando uma representação baseada em cadeias binárias, é claro que o espaço real coincide com o espaço conceptual. O espaço efectivo depende dos operadores. Se usarmos mutação por comutação de um bit e cruzamento de um ponto, todos os indivíduos se ligam entre si de modo directo ou, mais comum, de modo indirecto. O espaço efectivo é o mesmo que os dois anteriores.

(2) Passemos ao caso de procura do máximo (ou do mínimo) de uma função f: R^n --> R. O conjunto de todos os vectores de reais de tamanho será o nosso espaço conceptual. Admitamos que escolhemos uma representação directa, pelo que qualquer indivíduo é um vector de números reais de tamanho n. Podemos ser levados a concluir que, com esta representação, o espaço conceptual e o real coincidem. Mas tal não é verdadeiro, pois as limitações da representação de reais pelo computador fazem com que apenas uma parte possa ser representada. O espaço efectivo é também, em geral, um subconjunto do espaço real. Admitamos, por exemplo, que usamos apenas um operador de mutação que soma ou diminui uma unidade ao valor de cada componente do vector. Os vectores se podem ligar entre si, através deste operador, vão neste caso depender dos indivíduos da população inicial. É claro que por mutação de um elemento [1,2,3] nunca poderemos chegar a um elemento [1.5,2,3].

(3) Finalmente consideremos o problema das N-rainhas. O espaço conceptual é o espaço de todos os tabuleiros nXn com n rainhas colocadas, e tem por dimensão n ao quadrado levantado a n (maior que 4 milhares de milhões, para n = 8). Admitamos que decidimos representar as soluções considerando que existe apenas uma rainha em cada linha. Agora a dimensão será n levantado a n (maior que 16 milhões, para n=8). Se continuarmos com restrições, impondo que cada coluna só pode ter uma rainha, a dimensão do espaço de procura baixa (?) para “apenas” n factorial (igual a 40320, para n=8). Fixemos esta última representação. Fica claro que o espaço conceptual e o real são distintos, o primeiro contendo o segundo. E o espaço efectivo? Uma vez mais depende dos operadores usados. Admitamos escolher operadores, seja de mutação seja de cruzamento que garantem que os indivíduos resultantes da sua aplicação sejam permutações. Fica apenas por garantir que qualquer permutação pode ser obtida a partir de qualquer outra, ou seja, que independentemente da população inicial podemos chegar a qualquer ponto do espaço de procura real. Caso os operadores sejam liberais, isto é, possam produzir indivíduos que não sejam permutações, estamos a alterar o espaço de procura real, aumentando a sua dimensão. Temos que ter as consequências disso bem presente, em particular, determinar em que medida a qualidade destes indivíduos é medida de modo diverso do da qualidade dos indivíduos que são permutações.

Viáveis e Inviáveis

A questão que discutimos tem associada uma outra nuance, bem complexa. Em muitos problemas do mundo real, por exemplo problemas de escalonamento da produção, existem restrições que têm que ser respeitadas (e.g., uma máquina não pode estar a ser usada em dois processos diferentes). Se a nossa representação não incluir as restrições ao problema (no caso das n-rainhas que vimos acima, a representação por permutação permite incorporar duas das três restrições na representação) então é possível que o espaço real contenha aquilo a que chama soluções inviáveis. Nestes casos podemos optar por escolher os operadores de modo a que o espaço efectivo apenas contenha viáveis. No entanto, pode acontecer que o espaço das soluções viáveis esteja particionado, como se ilustra na figura abaixo.





Se por acaso todos os indivíduos pertencerem a um sub-espaço que não contém a solução óptima, e não for possível passar de um sub-espaço viável para outro igualmente viável, estaremos em maus lençóis. Muitas vezes permite-se então que se guardem soluções inviáveis, na esperança que estas nos ajudem a transitar entre espaços. Agora, como já foi referido mais acima, ficamos com a questão de saber como avaliar esses indivíduos. Uma técnica consiste em penalizá-los, por forma a que tenham baixa qualidade mas se possam manter temporariamente na população.

Ver para crer

Perceber como um algoritmo está a funcionar é muito importante. No caso dos algoritmos de que vamos tratar ainda mais. Visualizar o modo como a qualidade das soluções evolui ao longo do tempo ajuda imenso a perceber o que se está a passar. Para nos ajudar nisso vamos recorrer ao módulo matplotlib.

Vejamos como podemos fazer a visualização nos casos simples que nos interessam.

import matplotlib.pyplot as plt

def mostra_dados(dados):
    plt.title('')
    plt.xlabel('Generation')
    plt.ylabel('Cost')
    x = range(len(dados))
    y = dados
    plt.plot(x,y,'r-o')
    plt.show()

Se consultar o manual do matplotlib verá que pode fazer muito mais.

Agora só precisamos de alterar o programa principal, coligindo os valores ao longo das iterações para, no final, proceder à sua visualização.

def meu_padrao(problema,max_iter):
    candidato = prima_sol(problema)
    custo_candidato = qualidade(candidato)
    dados = [custo_candidato]
    for  i in range(max_iter):
        novo_candidato = proxima_solucao(problema, candidato)
        custo_novo_candidato = qualidade(novo_candidato)
        if melhor(custo_novo_candidato,custo_candidato):
            candidato = novo_candidato
            custo_candidato = custo_novo_candidato  
        dados.append(custo_candidato)
    mostra_dados(dados)
    return candidato,custo_candidato

Correndo o programa eis um resultado possível.




Como veremos mais adiante, porque se trata de algoritmos estocásticos, o normal é correr o algoritmo várias vezes, tipicamente um número maior ou igual a 30, e apresentar os valores médios. Eis um resultado possível.



É uma curva típica: uma melhoria rápida na fase inicial e depois maior dificuldade em nos aproximarmos do óptimo. A alteração ao código envolve o programa principal que é alterado por forma a devolver o resultado de cada execução e o modo como se faz a preparação dos dados. Vejamos como.

if __name__ =='__main__':
    resultados = []
    for i in range(30):
        padrao, dist, dados = meu_padrao((4,5),100)
        resultados.append(dados)
    
    valores = zip(*resultados)
    final = [ sum(gera)/float(len(gera)) for gera in valores]
    mostra_dados(final)

Notar o modo como se usa a função zip.

Recordando uma aula

Na aula extraordinária (?) da semana passada sugeri um pequeno exemplo para ilustrar como é simples a implementação de um algoritmo de optimização. O exemplo envolvia o reconhecimento de padrões simples, no caso caracteres; o algoritmo era um trepa colinas básico -- estocástico, uma solução de cada vez, busca local. O desenvolvimento foi feito em tempo real o que permitiu irmos discutindo alguns aspectos da implementação, incluindo os erros. A abordagem seguida para o desenvolvimento foi do topo para a base, usando camadas de abstracção.

Começámos então pelo programa principal.

def meu_padrao(problema,max_iter):
    candidato = prima_sol(problema)
    custo_candidato = qualidade(candidato)     
    for  i in range(max_iter):
        novo_candidato = proxima_solucao(problema, candidato)
        custo_novo_candidato = qualidade(novo_candidato)
        if melhor(custo_novo_candidato,custo_candidato):
            candidato = novo_candidato
            custo_candidato = custo_novo_candidato   
    return candidato,custo_candidato

Com esta aproximação o que abstraímos foram:
o problema
a geração aleatória de soluções
a geração dos vizinhos
a determinação da qualidade das soluções

Começando pelo problema. Resolvemos considerar uma discretização dos caracteres, ou seja, usar matriz de mXn células, cada uma podendo assumir o valor preto ou o valor branco. A representação correspondente em Python, mais simples, é uma lista de listas. Assim, por exemplo, o número um será representado pela lista.

[[0,0,1,0,0],[0,1,1,0,0],[0,0,1,0,0],[0,0,1,0,0]]

Definido o problema e a representação, podemos passar às restantes questões.

Geração aleatória de soluções

A estratégia é a de construir a matriz que represetna o caractere, linha a linha.

 def prima_sol_b(problema):
    mat = []
    for lin in range(problema[0]):
        linha = []
        for col in range(problema[1]):
            linha.append(random.randint(0,1))
        mat.append(linha)
    return mat

Geração de um vizinho

Escolhemos um índice de forma aleatória e mudamos o seu valor. Recorre-se a uma cópia da solução

def proxima_solucao(problema,solucao):
    linha = random.randint(0,problema[0]-1)
    coluna = random.randint(0,problema[1]-1)
    pixel = solucao[linha][coluna]
    novo_pixel = (pixel + 1) % 2
    nova_solucao = copy.deepcopy(solucao)
    nova_solucao[linha][coluna] = novo_pixel
    return nova_solucao

Note-se o recurso ao módulo copy e à operação deepcopy. Isto é necessário para evitar a modificação, por efeito lateral, da solução dada como referência, visto estarmos a usar listas de listas.

Qualidade

A qualidade de um indivíduo é dada pelo número de posições em que é diferente do caractere de referência.

def qualidade(solucao):
    padrao = [[0,0,1,0,0],[0,1,1,0,0],[0,0,1,0,0],[0,0,1,0,0]]
    return distancia(solucao,padrao)

def distancia(candi_1, candi_2):
    return sum([1 for linha in range(len(candi_1)) for coluna in range(len(candi_1[0]))\
                if candi_1[linha][coluna] != candi_2[linha][coluna]])

Acrescentamos a tudo isto a comparação de qualidade das soluções. Como calculamos a qualidade usando a distância é óbvio que quanto menor esse valor, melhor:

def melhor(quali_1, quali_2):
    return quali_1 < quali_2

E pronto. Já está. É só juntar tudo. Ou talvez não! Será que não estamos ainda satisfeitos? É um facto que as funções para geral soluções aleatoriamente e calcular a distância seguem um padrão que nos leva a pensar em listas por compreensão. Feita a substituição, aqui fica o código final.

import random
import copy

def meu_padrao(problema,max_iter):
    candidato = prima_sol(problema)
    custo_candidato = qualidade(candidato)     
    for  i in range(max_iter):
        novo_candidato = proxima_solucao(problema, candidato)
        custo_novo_candidato = qualidade(novo_candidato)
        if melhor(custo_novo_candidato,custo_candidato):
            candidato = novo_candidato
            custo_candidato = custo_novo_candidato   
    return candidato,custo_candidato


def prima_sol(problema):
    sol = [[random.randint(0,1) for col in range(problema[1])] for linha in range(problema[0])]
    return sol 

 
def proxima_solucao(problema,solucao):
    linha = random.randint(0,problema[0]-1)
    coluna = random.randint(0,problema[1]-1)
    pixel = solucao[linha][coluna]
    novo_pixel = (pixel + 1) % 2
    nova_solucao = copy.deepcopy(solucao)
    nova_solucao[linha][coluna] = novo_pixel
    return nova_solucao


def qualidade(solucao):
    padrao = [[0,0,1,0,0],[0,1,1,0,0],[0,0,1,0,0],[0,0,1,0,0]]
    return distancia(solucao,padrao)


def melhor(quali_1, quali_2):
    return quali_1 < quali_2

def distancia(candi_1, candi_2):
    return sum([1 for linha in range(len(candi_1)) for coluna in range(len(candi_1[0]))\
                if candi_1[linha][coluna] != candi_2[linha][coluna]])

                
def mostra(padrao):
    for linha in padrao:
        print linha, '\n'

if __name__ =='__main__':   
    padrao_ref = [[0,0,1,0,0],[0,1,1,0,0],[0,0,1,0,0],[0,0,1,0,0]]
    pad,dist = meu_padrao((4,5),50)
    mostra(padrao_ref)
    print 20 * '+'
    mostra(pad)
    print dist
    print qualidade(pad)

Retomar o caminho

Este blogue tem tido pouca actividade. Agora que começa uma nova edição da cadeira, vamos tentar manter uma intervenção mais regular. A ideia é a mesma: exemplos, desafios, ideias e outras coisas mais, tudo relacionado com a computação de inspiração biológica. Contribuições são bem vindas.