Números de João Brandão

Existem problemas de enunciado simples mas que se revelam bem complexos quando se trata de os resolver. Um exemplo disso é o problema conhecido pelo nome de Os números de João Brandão. Reza assim. Consideremos o conjunto dos números inteiros entre 1 e um dado n. Formar o seu sub-conjunto de maior cardinalidade em que nenhum dos elementos é igual à média da soma de outros dois membros presentes no sub-conjunto. Por exemplo, para n=5 uma solução possível é {1,2,4,5}. É claro quanto maior for o valor de n mais complexa é a questão.

Qual será a melhor forma de o resolver? Um método simples é o de enumerar todos os sub-conjuntos de inteiros que é possível construir e seleccionar os que não violam a restrição indicada. Depois é fácil: basta escolher os de maior cardinalidade. Mas sabemos quantos sub-conjuntos podem existir para um dado n, para determinar se este método de força bruta método é exequível? Se pensarmos um pouco chegamos à conclusão que esse número é igual à soma de todas as combinações de n tomadas i a i, com i a variar entre 1 e n. Feitas as contas esse valor é igual a 2 levantado a n. Um número muito grande para valores de n razoáveis...

Um alternativa a esta questão da dimensão é usar um método heurístico, por exemplo, um dos métodos estocásticos de estado único dados nas aulas. Não temos a certeza de encontrar uma solução óptima, mas pelo menos teremos um valor ... antes de morrermos. Dito isto, uma primeira questão que se coloca é a de como representar o problema e de como definir a qualidade de uma solução. A primeira resolve-se facilmente: um vector binário, em que um 1 na posição i significa que o número i faz parte da solução. Quanto à qualidade temos duas coisas a considerar: (1) não existirem números que violem a restrição enunciada e, (2) ter o tamanho maior possível. Na solução apresentada iremos considerar estes dois aspectos. Vejamos como.

Dada uma solução na forma de um vector binário a qualidade será tanto melhor quanto maior for o número de números presentes, e isso é dado pela quantidade de uns no vector. Por outro lado, temos que calcular a quantidade de violações da restrição. Daí o código:

def evaluate(indiv,comp):
    alfa = 1
    beta = 1.5
    return alfa * sum(indiv) - beta * viola(fenotipo(indiv),comp)

Esta solução faz intervir dois parâmetros, alfa e beta, que pesam do ponto de vista relativo os dois aspectos (tamanho e violações). Notar que quanto maior for o valor achado, melhor!

Para calcular o número de violações vamos usar o conjunto dos inteiros correspondente à cadeia binária (função fenótipo):

def fenotipo(indiv):
    fen = [i+1 for i in range(len(indiv)) if indiv[i] == 1]
    return fen

O algoritmo que calcula o número de violações parte de uma ideia simples: um número está a mais se os inteiros à mesma distância dele também estiverem no conjunto. Por exemplo, se eu tiver o número 4, então a presença simultânea do par (3,4), (2,6) ou (1,7) provoca uma violação cada. Usando esta ideia chegamos ao código:

def viola(indiv,comp):
    # count violations of constraints
    v = 0
    for elem in indiv:
     limite = min(elem-1,comp - elem)
     vi = 0
     for j in range(1,limite+1):
      if ((elem - j) in indiv) and ((elem+j) in indiv):
       vi += 1
     v += vi
    return v

Resolvido esta primeira parte da busca de uma solução, a seguir escolhemos um algoritmo, por exemplo, o Trepa Colinas em que o vizinho que consideramos para comparar com a solução candidata corrente é o melhor vizinho deste último:

def jb_hc(problem_size, max_iter,fitness):
    candidate = random_indiv(problem_size,problem_size)
    cost_candi = fitness(candidate)
    for i in range(max_iter):
        next_neighbor = best_neighbor(candidate,fitness)
        cost_next_neighbor = fitness(next_neighbor,problem_size)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor        
    return candidate

Notar como é feita a comparação, pois, como dissemos, quanto maior o valor dado pela função que calcula o mérito melhor. E como calculamos o melhor vizinho? em problemas discretos em que o número de vizinhos é finito e pequeno a resposta é fácil:

 def best_neighbor(individual, fitness):
    best = individual[:]
    best[0] = (best[0] + 1) % 2
    for pos in range(1,len(individual)):
        new_individual = individual[:]
        new_individual[pos]= (individual[pos] + 1) % 2
        if fitness(new_individual) > fitness(best):
            best = new_individual
    return best

Na solução apresentada vamos fabricando de modo ordenado os vizinhos de uma solução candidata, considerando que estes são aqueles que diferem em apenas uma posição dessa solução. à medida que geramos uma nova solução vamos actualizando, se for caso disso, o melhor vizinho. Juntando tudo, aqui fica o código completo. Agora é só experimentar para diferentes valores de n. Uma nota: quando n é igual a 100, se conseguir uma solução de tamanho 27 ... diga-me!

import random

# Basic Hill Climbing

def jb_hc(problem_size, max_iter,fitness):
    candidate = random_indiv(problem_size,problem_size)
    cost_candi = fitness(candidate)
    for i in range(max_iter):
        next_neighbor = best_neighbor(candidate,fitness)
        cost_next_neighbor = fitness(next_neighbor,problem_size)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor        
    return candidate
 
    
# Random Individual
def random_indiv(size):
    return [random.randint(0,1) for i in range(size)]

# Best neighbor
def best_neighbor(individual, fitness):
    best = individual[:]
    best[0] = (best[0] + 1) % 2
    for pos in range(1,len(individual)):
        new_individual = individual[:]
        new_individual[pos]= (individual[pos] + 1) % 2
        if fitness(new_individual) > fitness(best):
            best = new_individual
    return best


# Fitness for JB
def evaluate(indiv,comp):
    alfa = 1
    beta = 1.5
    return alfa * sum(indiv) - beta * viola(fenotipo(indiv),comp)

def viola(indiv,comp):
    # count violations of constraints
    v = 0
    for elem in indiv:
     limite = min(elem-1,comp - elem)
     vi = 0
     for j in range(1,limite+1):
      if ((elem - j) in indiv) and ((elem+j) in indiv):
       vi += 1
     v += vi
    return v
    
# Auxiliar

def fenotipo(indiv):
    fen = [i+1 for i in range(len(indiv)) if indiv[i] == 1]
    return fen

    
if __name__ == '__main__':
 # For test purposes: beware of the time it may takes...
    res = fenotipo(jb_hc(100,10000,evaluate))
    quali = viola(res, len(res))
    print(quali,'\n', len(res), '\n',res)