Programming Pearls: Parte 1

Janeiro 28, 2011

Estou lendo o livro “Programming Pearls” do Jon Bentley. O livro consiste de um agregado de diversas colunas (publicadas no Communications of ACM) abordando aspectos práticos de programação. Uma curiosidade: Jon Bentley é autor do conhecido algoritmo de geometria computacional para encontrar as interseções de um conjunto de segmentos de reta, o algoritmo de Bentley-Ottmann.

O livro é dividido em 15 colunas, agrupadas em 3 partes. Pretendo fazer posts comentando as partes que mais gostei em cada parte. Começo pela parte I, que abrange as colunas 1 a 5.

Coluna 1: Cracking the Oyster

Nesse capítulo o autor apresenta um problema de ordenar um arquivo com entradas com valor de 1 a únicas, com limitação de memória.

Podemos usar um algoritmo de ordenação linear. A sacada é usar um vetor de bits. O i-ésimo bit é 1 se o valor i está no arquivo. Logo são necessários bits ao invés dos se fôssemos usar inteiros de 32-bits.

Coluna 2: Aha Algorithms

Aqui ele apresenta alguns algoritmo de “sacada”. O que achei mais interessante foi o seguinte: rotacionar um vetor de tamanho n, i posições à equerda, em tempo O(n) e espaço adicional O(1).

A ideia que eu tive inicialmente satisfaz as restrições de tempo-espaço, mas é bem complicada. Ao rotacionar k posições, sabemos que v[0] receberá v[i] e v[i] receberá v[2i] (sendo que talvez 2i eventualmente dê a volta do vetor). Seguimos trocando v[ki] com v[(k+1)i] até que alguma hora voltamos no 0.

O problema é que isso pode não percorrer todos os elementos. Considere um vetor v com 6 elementos e que queremos rotacioná-lo em 2 posições.

Aqui temos v[0] <- v[2] <- v[4] <- v[0]. Precisamos fazer também iniciando em v[1], ou seja, v[1] <- v[3] <- v[5] <- v[1]. A pergunta é: como saber quando um elemento não foi ainda coberto, se só podemos usar um número constante de espaço extra?

Para isso, temos a seguinte proposição: Seja m = mdc(i, n). Então um ciclo iniciado no elemento v[x] cobre a posição j se e só se

Esboço da prova: Se j é visitado, então para algum k inteiro não-negativo, ou seja, x + ki = qn + j, para algum q inteiro não-negativo. Como por definição m divide i e n, o resto de j por m dá x. Se j não é visitado, dá para mostrar que o resto de j por m é diferente de x.

No exemplo acima, mdc(i=2, n=6) = 2. Logo, o ciclo iniciado em 0 vai cobrir 0, 2 e 4 e o iniciado em 1 vai cobrir 1, 3 e 5.

Um corolário dessa proposição é que podemos percorrer m ciclos, começando de 0, 1, …, m-1, que todos os elementos serão contemplados.

O autor apresenta uma solução extremamente simples e elegante. Primeiro, divida o vetor V em 2 partes. De 0 a i-1 e de i a n-1. Chame-as respectivamente de A e B.


Vetor original

Rotacionar o vetor correspondente à figura acima resulta na figura abaixo.

Vetor rotacionado

Que é simplesmente trocar os blocos A e B de lugar. Para fazer isso, considere a operação que inverte a ordem dos elementos em um vetor. É fácil fazer in-place e em O(n).

Vamos denotar por A’ e B’ os vetores A e B invertidos. Ao inverter o vetor V = AB, teremos B’A’ . Ora, para chegar em BA, basta inverter separadamente A’ e B’!

Coluna 3: Data Structures Program

Esse capítulo discute o uso de estruturas de dados, principalmente vetores, não para melhorar a eficiência dos algoritmos, mas para diminuir o tamanho do código. Não é apresentada nenhum técnica interessante, é mais uma dica mesmo.

Coluna 4: Writing Correct Programs

Aqui o autor faz uso de assertivas (os assert‘s em C e python), incluindo pré e pós condições para checar invariantes de laços e funções. Para isso ele usa como exemplo o algoritmo de busca binária, que é relativamente simples, mas muito tricky de implementar. (Tão tricky que o caderno de algoritmos do time do ITA para o mundial da maratona de 2008 tinha umas 4 implementações de busca binária).

Discutindo sobre quais assertivas usar para ter certeza de que a busca binária termina, ele apresenta dois exercícios interessantes. Primeiro, mostrar que o seguinte algoritmo termina:

while x != 1:
    if x é par:
        x = x/2
    else:
        x = 3x + 1

Esse exercício é praticamente uma pegadinha, pois a resposta do exercício é conhecida como a conjectura de Collatz, um problema em aberto.

O outro é o seguinte: Dado um saco com W feijões brancos e B feijões pretos, repita o seguinte procedimento: escolha dois feijões aleatórios desse saco. Se tiverem a mesma cor, joga fora ambos os feijões e devolve um feijão preto ao saco. Se forem de cor diferente, devolve um feijão branco. Pede-se para mostrar que isso eventualmente terminar e caracterizar o resultado em função de W e B.

Coluna 5: A small matter of programming

Esse capítulo parece ser um complemento do anterior. Além de assertivas, o autor introduz testes automáticos e medidas de tempo ao código da busca binária.

Ele chama esses componentes adicionais de scaffolding, que faz alusão aos andaimes que se usam em construção de prédios e depois são removidos. Ao contrário da construção civil, porém, pode ser uma boa prática deixar esses componentes em código de produção (desde que o overhead não seja muito grande).

Anúncios

Branch and Bound: Seleção de variável

Janeiro 21, 2011

Um aspecto que pode influenciar o desempenho de um algoritmo de Branch and Bound é a escolha da variável que sofrerá branching. Para relembrar, em cada nó da árvore de exploração do Branch and Bound, resolve-se a relaxação linear do modelo, o que, na maioria das vezes, retorna uma solução fracionária.

A ideia é então escolher uma variável com valor fracionário por exemplo, e criar dois novos modelos, um com a restrição e outro com a restrição .

A questão é, qual a melhor variável para se escolher? Uma regra bem simples é escolher aquela com valor “mais fracionário”, ou seja, cuja parte decimal seja mais próxima de 0.5, desempatando arbitrariamente. Porém há pelo menos duas técnicas mais sofisticadas que iremos apresentar a seguir.

Strong Branching

O strong branching é basicamente o seguinte. Para cada variável, calculamos o valor da relaxação linear do PL com a restrição e do PL com a restrição , chamando-os respectivamente de e . Observe que ambos são limitantes duais para o problema. Calculamos uma soma ponderada desses dois valores e escolhemos a variável com valor mais apertado (ou seja, menor valor se o proplema é de maximização ou maior se é de minimização).

O problema dessa abordagem é o custo. Resolver duas relaxações lineares para cada variável em cada nó é bastante custoso. Uma abordagem aproximada consiste em pegar apenas as N variáveis “mais fracionárias” para fazer tal cálculo. Além do mais, executar no máximo K iterações do algoritmo dual-simplex e não até encontrar a solução ótima.

Pseudo-custos

Outra alternativa mais barata é o uso de pseudo-custos, que são uma estimativa dos valores e do strong branching. Essa estimativa é essencialmente a média obtida dos valores das relaxações lineares ao longo da execução do algoritmo ao fazer o branching das variáveis. Assim, se em um dado momento o algoritmo escolhe a variávei i para sofrer branching, analisamos os valores das relaxações lineares de seus filhos e atualizamos essas médias para i.

Para termos um valor inicial das estimativas de e , podemos aplicar o strong branching apenas na raíz.

Prioridades

Há vezes em que sabemos, por teoria, que fazer branching em uma certa classe de variável primeiro, tende a gerar árvores de busca de menor altura. Logo, um vetor de prioridades de cada variável fornecido como parte do modelo pode ser usado para guiar a escolha do algoritmo de Branch and Bound.

Como veremos a seguir, a biblioteca XPRESS trabalha com o conceito de prioridades e pseudo custos.

Seleção da variável no XPRESS

Segundo o manual, a biblioteca XPRESS trabalha com dois parâmetros para cada variável: a prioridade (o padrão é 500) e os pseudo-custos (up e down). O trecho a seguir, extraído desse manual, explica o algoritmo:

Each global entry has a priority for branching, either the default value of 500 or one set by the user in the directives file. A low priority value means that the variable is more likely to be selected for branching. Up and down pseudo costs for each global entity can be specified, which are estimates of the per unit degradation of forcing the entity away from its LP value.

The Optimizer selects the branching entity from among those entities of the most important priority class which remain unsatisfied. Of these, it takes the one with the highest estimated cost of being satisfied (degradation).

O algoritmo divide as variáveis em classe de prioridades. Pertencem à mesma classe aquelas variáveis que possuem mesma prioridade. O algoritmo então escolhe a classe com menor valor de prioridade que tenha pelo menos uma varíavel fracionária. Para escolher a variável dentro dessa classe, ele utiliza os pseudo-custos como critério de desempate.

Ele define o valor de degradação como uma combinação dos pseudo-custos up e down de cada variável. Qual a conta feita depende do parâmetro VARSELECTION, mas algumas delas são: min(up, down), up + down e min(up, down) + max(up, down). Então, ele seleciona a variável com maior valor de degradação.

Referências

[1] Notas de aula sobre Programação Inteira por Jeff Linderoth

[2] Encyclopedia of Optimization, 2nd Edition – Chapter: Integer Programming: Branch and Bound Methods, Branching Variable Selection


Topcoder: SRM 492

Janeiro 14, 2011

Topcoder é uma empresa que realiza diversos tipos de competições de programação. A mais famosa, de longe, é o SRM (Single Round Match), uma competição rápida de algoritmos. Além dessa, há competições de design e desenvolvimento de software, além das competições longas de algoritmos (Marathon Matches). Veja o site para mais informações.

O SRM é dividido em 3 fases: código, desafio e testes do sistema. A fase de código dura 75 minutos e consiste em resolver 3 problemas estilo maratona de programação. A diferença é que você submete no escuro, ou seja, sem saber se seu programa passou nos casos de testes do sistema. A fase de desafio te dá uma oportunidade de ganhar uns pontos extras fornecendo casos de teste que possivelmente quebre o código do oponente e dura 10 minutos. Depois disso, os programas de cada competidor são submetidos aos testes do sistema. Tem vários detalhes sobre pontuação que podem ser encontrados aqui .

Eu participo dos SRM’s desde 2007, mas sempre de maneira esporádica. Agora que me aposentei da maratona, decidi continuar praticando para não enferrujar (embora eu ainda programe bastante no meu mestrado). Uma coisa que eu deveria fazer mas nunca tive disciplina suficiente, é o post mortem da competição, ou seja, estudar os editoriais e aprender a fazer os problemas que não consegui durante a competição.

Minha ideia é tentar usar o blog para escrever sobre tais problemas. Sem mais delongas, vou apresentar uma breve descrição e a solução dos problemas que não consegui resolver.

Time Traveling Tour

Dado um grafo ponderado e não-orientado G e uma lista de cidades L, deseja-se encontrar um passeio iniciado no vértice 0, que visite as cidades na ordem em que aparecem em L (o passeio pode conter cidades intermediárias que não necessariamente são visitadas). Há um diferencial porém: nesse passeio você pode viajar no tempo e voltar a uma cidade pela qual já passou, com custo zero. Queremos o passeio de menor custo. Link para o problema.

Se não houvesse essa viagem do tempo o problema seria fácil: ache o menor caminho entre todos os pares de vértices e some a distância entre cidades consecutivas de L.

A solução do editorial, descreve um algoritmo de programação dinâmica. O estado é definido por uma tripla (C, i, j), que representa o menor custo de se sair da cidade C e visitar as cidades L[i], …, L[j]. A resposta do problema é dada por (0, 0, |L|-1).

A sacada dessa abordagem é que não importa qual o caminho feito até chegarmos em C. Vamos inverter o pensamento: podemos ir a uma cidade D e então voltar para a cidade C com custo zero. Logo, a partir de uma dada cidade C, iremos usar outras cidades para visitar trechos de i, … j.

Considere o exemplo abaixo, onde as cidades são identificadas por letras e os números indicam a ordem em que devem ser visitadas. A cidade de partida é ‘A’.

Exemplo: letras são rótulos do vértices; números sobre os vértices representam a ordem em que devem ser visitados; os números em azul sobre as arestas representam os custos.

Podemos ver que uma solução ótima é a seguinte: vai de A a B com custo 3 e visita B, vai então até D com custo 2. Volta no tempo para A e visita C com custo 3. Volta no tempo para A e vai até E com custo 8. O total fica 16.

Em termos da tripla, poderíamos representar (A, 0, 3) como {dist(A, B) + (B, 0, 1)} + {dist(A, C) + (C, 2, 2)} + {dist(A, B) + (B, 3, 3)}. Por sua vez, (B, 0, 1) = (B, 0, 0) + {dist(B, D) + (D, 1, 1)}. A base da recursão é quando i = j ou i > j. No primeiro caso, (C, i, j) = dist(C, L[i]). No segundo significa que já visitamos todas as cidades então (C, i, j) = 0.

Uma observação importante é que {dist(A, B) + (B, 0, 1)} + {dist(A, C) + (C, 2, 2)} + {dist(A, B) + (B, 3, 3)} pode ser reescrita como

(1) {dist(A, B) + (B, 0, 1)} + (A, 2, 3),

pois (A, 2, 3) será recursivamente decomposto até que eventualmente a primeira forma apareça.

Isso reduz nosso trabalho a encontrar uma cidade (no exemplo é B) e um índice k, tal que essa cidade vai visitar as cidades i, …, k. A recorrência então fica:

(C, i, j) = INF
Para cada cidade D :
  Para cada índice k = i ... j :
      (C, i, j) = min {(C, i, j), 
                    dist(C, D) + (D, i, k) + (C, k+1, j)

Há dois problemas: essa recorrência pode ser muito ineficiente para os limites do problema, já que o número de operações fica O(50^5). Outro problema é quando k = j. Nesse caso teremos uma recorrência cíclica, pois chamaremos (D, i, j) que por sua vez pode chamar (C, i, j) novamente.

A solução do competidor those trata esses dois problemas. Note que só precisamos considerar o caso k = j porque pode existir uma cidade para a qual não se volta com a máquina do tempo. Então definimos uma segunda tripla {C, i, j} que corresponde ao menor custo de se visitar as cidades i e j, só que necessariamente voltando no tempo para C pelo menos uma vez.

Como sabemos que {C, i, j} voltará ao menos uma vez usando a máquina do tempo, podemos quebrá-lo em duas partes (antes e depois dessa volta):

(2) {C, i, j} = min{(C, i, k) + (C, k + 1, j)} para algum k = i … j – 1

Agora (C, i, k) e (C, k + 1, j) podem ou “delegar” a tarefa de visitar as cidades i a j para outras cidades ou ele mesmo quebrar em mais pedaços (que seria o equivalente a delegar a tarefa para a própria cidade C).

A recorrência de (C, i, j) fica assim:

(3) (C, i, j) = min{ {D, i, j} + dist(C, D)} para alguma cidade D

Observe que essa abordagem elimina a recorrência cíclica e de quebra reduz a complexidade para O(50^4). O trecho que calcula essa recorrência é o seguinte:

i64 eval2(int C, int i, int j){ 
  if (pd2[C][i][j] != -1) return pd2[C][i][j]; 
  pd2[C][i][j] = INFLL; 
  for(int k=i; k<j; k++)
    pd2[C][i][j] = min(pd2[C][i][j], eval(C, i, k)+ 
                        eval(C,k+1, j));
  return pd2[C][i][j]; 
} 
i64 eval(int C, int i, int j){ 
  if (pd[C][i][j]!=-1) return pd[C][i][j]; 
  if (i == j) return dist[C][L[i]]; 
  pd[C][i][j] = INFINITY; 
  for(int D=0; D<n; D++) 
    pd[C][i][j] = min(pd[C][i][j], dist[C][D]+
                      eval2(D, i, j)); 
  return pd[C][i][j]; 
} 

Sendo que a resposta do problema é dada por eval(0, 0, |L|-1).

A minha ideia era postar uma explicação para o problema mais difícil também, o TimeTravellingGogo. Porém, achei tão difícil e levei tanto tempo para pegar a ideia do problema TimeTravellingTour que desanimei. O próprio autor dos problemas disse que eles estavam bem difíceis. Tanto que somente uma pessoa acertou o problema de 1000 pontos.

É incrível que 31 pessoas tenham acertado o problema TimeTravellingTour! Isso quer dizer que eles tiveram a ideia e implementaram corretamente em 75 minutos! Eu lendo o editorial levei umas 2 semanas só para entender a ideia :/


Biblioteca OpenMP

Janeiro 7, 2011

Depois de estudar um pouco a biblioteca MPI e aprender suas principais funções, decidi investir um tempo na biblioteca OpenMP. Como eu já tinha comentado em outro post, a diferença principal entre as duas é que o MPI usa o paradigma de memória distribuída enquanto no OpenMP a memória é compartilhada.

Dispor de memória compartilhada pode facilitar bastante o trabalho. Um exemplo é a divisão balanceada das iterações entre os processadores e com trabalho totalmente distribuído (no referido post a distribuição era feita através de um processador mestre). Com o MPI chegamos a uma solução com comunicação não-bloqueante, mas que executava um número de iterações maior do que o previsto.

Com o OpenMP podemos manter um pool central de iterações. Conforme cada processador termine de executar seu bloco de iterações, ele adquire exclusividade sobre esse pool, atualiza o número de iterações restantes e libera o recurso.

Fiz a simulação de uma iteração usando a função sleep com tempos aleatórios.

Geração de números aleatórios

Em geral, um requerimento para reproducibilidade de experimentos em que se usam números aleatórios é que toda vez que executamos o código, a sequência gerada seja a mesma. Para isso a biblioteca padrão stdlib oferece a função srand, que recebe como parâmetro uma semente. Se a semente for fixada, a sequência terá essa propriedade.

Quando estamos trabalhando com paralelização de código, uma propriedade desejável é que a união das sequências geradas em paralelo correspondam à sequência que seria gerada sequencialmente. Por exemplo, se com um processo geramos 6 números aleatórios entre 1 e 10, {6, 6, 8, 9, 7, 9}, então essa mesma sequência deveria ser gerada para 2 processos. Poderia ser algo como {6, 8, 7} e {6, 9, 9} (na verdade depende da ordem em que cada processador invocou o gerador – o importante é que a união das sequências seja igual à sequência original).

Com isso, os resultados obtidos serão essencialmente os mesmos entre a versão sequencial e a versão paralela, restando comparar apenas o tempo de execução.

O problema do gerador da biblioteca padrão C é que não temos acesso ao estado do gerador. Assim, não dá para garantir que quando cada processador gerar um número aleatório, o gerador esteja com o mesmo estado em todos os processadores. A biblioteca Boost possui uma implementação de geradores de números aleatórios, com a vantagem de que o gerador é um objeto ao qual temos acesso.

Assim, podemos colocá-lo como variável compartilhada entre os processadores e colocar a geração do número aleatório dentro de uma região crítica do Open MP, que é uma região onde apenas um processador pode executar por vez (quem tentar acessar fica bloqueado até que o processador na região crítica termine). Com isso garantimos que o estado do gerador está sempre atualizado para todos os processadores. Essa abordagem pode ser conferida na seção de código abaixo.

Resultados computacionais

Rodei o código abaixo em uma máquina quad-core e obtive os seguintes tempos de execução (medidos com a função time do linux), em função do número de threads utilizadas:

Tempos de execução

Observamos que a aceleração é quase linear. A terceira coluna da tabela mostra um tempo “normalizado” crescente com o aumento do número de processadores, mas isso se deve ao trecho não paralelizado do código. De forma geral, a economia de tempo foi bastante satisfatória. Porém, trata-se de um exemplo muito simples. Pretendo ainda aplicar essa estratégia em uma implementação mais complexa para ver os ganhos obtidos.

Código

#include <stdio.h>
#include <iostream>
#include <algorithm>

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>
#include <omp.h>

/* A otenção do número aleatório deve ser protegido */
int random(boost::mt19937 &gen, int lb, int ub){
    int val;
    #pragma omp critical
    {
        boost::uniform_int< > dist(lb, ub);
        boost::variate_generator<boost::mt19937&, 
                                 boost::uniform_int<> > r(gen, dist);
        val = r();
    }
    return val;
}
int iteration(boost::mt19937 &gen){
    int rtime = random(gen, 1, 10);
    sleep(rtime);
}
/* A atualização das iterações restantes deve ser 
 * uma area protegida */
bool get_iteration_block(int *iter){
    bool stat;
    #pragma omp critical
    {
        if (*iter <= 0){
            stat = false;
        }
        else {
            *iter = *iter - 1;
            stat = true;
        }
    }
    return stat;
}
int main (){

    boost::mt19937 gen;
    gen.seed(0);

    int iter = 20;  // numero inicial de iteracoes
    #pragma omp parallel shared(iter, gen) num_threads(2) 
    {
        while (get_iteration_block(&iter))
            iteration(gen);         
    }
    return 0;
}

Para compilar código OpenMP no linux podemos usar o gcc com a flag -fopenmp

g++ -fopenmp <arquivo.cpp>

No código acima, o gerador mt19937 é baseado no algoritmo gerador de números pseudo-aleatórios Mersenne Twister.