As Figurinhas dos Vingadores são Honestas? (parte 2)

No último post do blog, eu comprei vários pacotinhos de figurinhas dos Vingadores, para tentar responder à pergunta: as figurinhas são honestas? Após fazer dois testes de hipótese, eu concluí que as figurinhas possuem distribuição uniforme, então parece que sim. Mas tem uma pegadinha: para serem honestas, garantir que são uniformes não é o suficiente, ainda tem outro teste que precisamos fazer!

A Distribuição dos Pacotinhos

Para fazer a análise das figurinhas, eu comprei 106 pacotinhos na banca. Mas enquanto eu abria os pacotinhos, notei uma propriedade curiosa: todos os pacotinhos tinham sempre 4 figurinhas diferentes. Mas olha só, se os pacotinhos não possuem repetidas entre si, então as figurinhas podem até ser uniformes, mas não são independentes! Cada pacotinho tem quatro figurinhas, e se eu tirei a figurinha 42 em um pacotinho, então eu sei que as outras três não são a figurinha 42.

Mas pode ter sido só o acaso né? Vai ver eu dei sorte de comprar exatamente 106 pacotinhos que não tinham repetidas. Mas felizmente nós sabemos como testar se foi sorte ou não, basta fazer o teste do ##\chi^2##, que vimos no último post!

Acho que o jeito mais fácil de testar é contar quantas figurinhas diferentes tem dentro do pacotinho. Quantos pacotinhos possíveis existem com 4 figurinhas diferentes? Esse é fácil, o álbum tem 180 figurinhas, então são ##180\choose 4## pacotinhos. Com 3 diferentes precisa de mais cuidado, existem ##180\choose 3## combinações de figurinhas, mas você precisa contar qual dessas três diferentes está repetida (os casos são AABC, ABBC, ABCC). Então para 3 diferentes são ##3 {180\choose 4}## pacotinhos.

Com 2 diferentes também tem três casos: AAAB, AABB, ABBB, logo temos ##3{180\choose 2}## pacotinhos. Para uma única figurinha diferente, são apenas ##180\choose 1## pacotinhos. Agora podemos montar a tabela com quantos pacotinhos eu deveria ter visto, e quantos pacotinhos eu vi:

	1 diferente	2 diferentes	3 diferentes	4 diferentes
Total de pacotinhos	##180\choose 1##	##3{180\choose 2}##	##3{180\choose 3}##	##180\choose 4##
Total (numérico)	180	48330	2867580	42296805
Total (normalizado)	3.98e-6	1.07e-3	0.634	0.936
Esperado (total*106)	422e-6	0.113	6.72	99.163
Observado	0	0	0	106

Opa, nas duas primeiras categorias ##np\lt 5##, então eu não tenho amostras suficientes para usar o ##\chi^2##! Um jeito de consertar isso sem comprar mais pacotinhos é juntando as três primeiras categorias em uma só:

	Menos de 4 diferentes	4 diferentes
Esperado	6.837	99.163
Observado	0	106

Com isso conseguimos calcular o ##\chi^2## pela definição, chegando no valor de ##\chi^2=7.308##. Qual é o limite para rejeitar a hipótese? Se estivermos trabalhando com ##\alpha=5\%##, e um grau de liberdade, então o limite é:

InverseCDF[ChiSquareDistribution[1], 1 - 0.05]
>> 3.84146

Como ##7.308\gt 3.84146##, então temos que rejeitar a hipótese, e concluir que as figurinhas dos Vingadores não são independentes, e, como consequência, não são honestas, apesar de serem uniformes!

Felizmente para o colecionador, elas são desonestas mas estão jogando do seu lado. Se os pacotinhos diminuem a chance de ter repetidas, então o total esperado de pacotinhos que você precisa comprar é menor que no caso puramente honesto.

Dá para calcular o quanto você economiza com as figurinhas desonestas? Sure, vou mostrar dois métodos para fazer isso:

Cadeias de Markov

Esse primeiro método é mais complicado, mas vou começar por ele porque é mais genérico, e se aplica a mais casos. Imagine que eu já tenho 50 figurinhas no meu álbum, e eu abro mais um pacotinho. Se eu der azar, vão sair quatro figurinhas que eu já tenho. Se eu der sorte, pode sair até quatro que eu ainda não tinha. 

O método mais natural para modelar esse problema é com uma Cadeia de Markov: você considera o estado como sendo quantas figurinhas você tem no álbum, e as transições possíveis são quantas figurinhas você ainda não tinha no pacotinho.

Para montar a matriz de transição, nós precisamos da probabilidade de ir de um estado ##n## até um estado ##m##. Como cada pacotinho tem 4 figurinhas diferentes, então o ##m## precisa ser necessariamente ##n##, ##n+1##, ..., ##n+4##, e outras transições não são possíveis. 

Qual a chance de ir do estado ##n## para ##n+k##? Nós sabemos que as quatro figurinhas do pacotinho são diferentes entre si, então temos ##k## figurinhas novas e ##4-k## figurinhas que já tínhamos. A quantidade de pacotinhos que fazem essa transição é igual a ##{180-n\choose k}{n\choose 4-k}##, já que precisamos escolher ##k## figurinhas dentre que as não temos, e ##4-k## dentre as que temos. O total de pacotinhos é ##180\choose 4##, então a probabilidade total do estado ##n## para ##n+k## é ##{180-n\choose k}{n\choose 4-k}/{180\choose 4}## . Com isso podemos montar a matriz de transição ##M##:

$$ M= \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & 0.02 & 0.97 & 0 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & \approx0 & 0.04 & 0.95 & 0 & \dots & 0 &  0 & 0 \\ 0 & 0 & 0 & 0 & \approx0 & \approx0 & 0.06 & 0.93 & \dots & 0 &  0 & 0 \\ \vdots &\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0.95 &  0.04 & \approx0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 &  0.97 & 0.02 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 &  0 & 1 \\ \end{bmatrix}$$

Pela matriz podemos conferir a intuição: se você não tem nenhuma figurinha, o primeiro pacotinho vai garantidamente te dar 4 figurinhas novas, então o estado 0->4 na matriz vale 1, e todos os outros nessa linha são zero. Por outro lado, se está faltando só uma figurinha para completar, comprar um pacotinho novo dá um chance de 2% de conseguir a faltante. E, claro, se você tem todas as figurinhas, então a chance de você não conseguir mais nenhuma é 100%.

O que queremos calcular é quanto economizamos devido às figurinhas desonestas. Para isso, vamos calcular qual o valor esperado do número de pacotinhos que precisamos comprar para completar o álbum. Tendo a matriz de transição em mãos, isso é equivalente ao número médio de transições para ir do estado 0 (álbum vazio) ao estado 180 (álbum completo). E como calcular isso? Usando um algoritmo secreto!

Really, um algoritmo secreto. Não sei o motivo, mas você não acha esse algoritmo em lugar nenhum; não tem na wikipedia e nem nos livros de processos estocásticos que eu tenho aqui em casa. Então aproveitem para anotar porque é bem útil:

1. Dada a matriz de transição ##M##, calcule a matriz reduzida ##R##, retirando da matriz original todas as linhas e colunas associadas aos estados finais. (No nosso caso só a última linha e coluna são finais).
2. Calcule a matriz fundamental ##N## do processo, através da fórmula ##N=(I-R)^{-1}##.
3. O valor esperado do número de transições é a soma dos elementos da primeira linha da matriz ##N##.

Implementando o algoritmo acima com um scriptinho rápido no Mathematica, concluímos que o valor esperado do número de pacotinhos que você precisa comprar é ##257.98##. Como cada pacotinho custa R$ 1.20, então você precisa de R$309.58 em média para completar o álbum.

E quanto você gastaria com figurinhas 100% honestas? Esse cálculo eu já fiz aqui no blog anteriormente, se as figurinhas são uniformes e independentes então esse é o Coupon Collector Problem: para ##n## figurinhas a solução é ##n H(n)##, onde ##H(n)## é o número harmônico de ##n##. Essa fórmula dá o total de figurinhas, como elas são independentes, para o total de pacotinhos é só dividir por 4, chegando no total de ##259.78## pacotinhos, ou R$ 311.74, uma economia de 2 reais! (Eu falei que economizava, só não falei que era muito).

Embora o algoritmo secreto seja bacana e funcione em vários casos, ele é meio lento. Você precisa inverter uma matriz 181x181, o que não é muito fácil. O fato dela ser triangular e esparsa ajuda, mas tem um método melhor, que calcula o valor em O(n)!

Linearidade da Expectativa

Qual é o teorema mais surpreendente que você conhece? Só em probabilidade tem vários que são difíceis de acreditar à primeira vista, como o Monty Hall Problem e o Paradoxo do Aniversário. Mas um que sempre me espanta é a linearidade da expectativa (o valor esperado da soma é a soma dos valores esperados):

$$E[a+b]=E[a]+E[b]$$

Para mim a parte curiosa é que ele funciona mesmo se as variáveis ##a## e ##b## não forem independentes! Para o nosso problema isso ajuda bastante: o valor esperado para conseguir 180 figurinhas claramente é correlacionado com o valor esperado para conseguir 179 figurinhas, mas eu posso usar relações lineares entre eles sem medo.

Para essa análise eu acho mais fácil pensar invertido. Se estiver faltando 0 figurinhas para completar o album, quantos pacotinhos eu preciso? Nenhum né, ele já está completo. E se estiver faltando só uma figurinha? Aí depende do que sair no pacotinho, se sair a figurinha faltante, um pacotinho basta, se não sair, então você adiciona um pacotinho no seu valor esperado e repete o procedimento:

$$\begin{align*}\\E[1]&=p_0(1+E[0]) + p_1(1+E[1])\\ (1-p_1)E[1] &= p_1 + p_0(1+E[0]) \\ E[1] &=\frac{1}{1-p_1}\left(p_1+p_0(1+E[0])\right)\end{align*}$$

Olha só, eu sei que ##E[0]=0##; para calcular ##E[1]## eu preciso só do ##E[0]##; para o ##E[2]## eu preciso do ##E[1]## e ##E[0]##; e assim por diante. Como cada pacotinho tem no máximo 4 figurinhas que eu não tenho, então eu só preciso olhar quatro valores para trás na recorrência. Portanto, esse algoritmo é linear!

Preenchendo as probabilidades corretas, o algoritmo completo fica assim:

$$\begin{align*} \\ E[0] &= 0 \\ E[n] &= \frac{1}{1-{{180-n}\choose 4}/{180\choose 4}}\left(\frac{{180-n}\choose 4}{180\choose 4}+\sum_{1\le k \le 4}\frac{{n\choose k}{{180-n}\choose{4-k}}}{180\choose 4}E[n-k]\right) \\ \end{align*} $$

Essa implementação realmente é bem mais rápida que o anterior. A versão com matrizes no Mathematica leva 1s para rodar, enquanto que um scriptinho python com o algoritmo linear leva só 55ms!

Agora que temos os pacotinhos corretamente modelados, finalmente vamos conseguir bolar uma estratégia para minimizar o custo de completar o álbum. Mas vai ficar para o próximo post porque esse já ficou comprido demais :)

A Busca Por Bozo Binário

Se você trabalhar por tempo o suficiente com Computação, eventualmente uma das suas atribuições vai ser entrevistar candidatos para a sua equipe. Entrevistar é uma arte: existem inúmeros truques e macetes para descobrir se é esse sujeito na sua frente é mesmo o cara que você quer como seu parceiro. O Joel tem uma boa introdução sobre o tema que vale a pena ler.

Uma das habilidades que você precisa desenvolver rápido é aprender que o candidato não pensa como você. Em especial, quando ele responde uma pergunta de um jeito diferente daquele que você esperava, não quer dizer necessariamente que ele esteja errado. E esse caso eu já vi de perto, com o algoritmo que eu apelidei de Busca Por Bozo Binário.

Essa história aconteceu faz algum tempo. Eu estava na cantina conversando sobre entrevistas, quando alguém comentou que rejeitou um candidato porque ele não conseguiu implementar uma busca binária. Poxa, isso é grave mesmo! E qual o algoritmo errado que ele fez? A resposta foi essa:

"Sorteie um número de 1 até o tamanho do vetor, e verifique o elemento que tinha esse índice. Se for igual você achou, então retorna. Se for menor, repete para o lado esquerdo; se for maior, repete para o lado direito."

Isso é errado, me disseram, porque o pior caso é ##O(n)##. De fato, se o elemento que você quer achar é o primeiro, e o gerador de aleatórios retornar sempre o último, ele vai passar pelo vetor inteiro.

Mas nesse ponto eu tive que interromper: "Peraí! Você explicitou que queria ##O(\log n)## no pior caso? Se você não falou nada no enunciado, ele pode ter entendido que era ##O(\log n)## no caso médio, e aí o algoritmo dele está certo!"

Na hora ninguém botou fé que esse algoritmo de Busca por Bozo Binário era de fato ##O(\log n)##, então eu tive que voltar para casa e escrever a demonstração. Se você tem medo de Matemática Discreta, pule a caixa azul:

Para a Busca por Bozo Binário, nós queremos calcular qual é o valor médio do número de comparações que ele realiza. E como você começa uma análise de caso médio?

Primeiro você define as características da entrada: eu vou assumir que todos os casos estão distribuídos uniformemente. Depois, é só usar a fórmula do valor esperado: você soma o número de comparações em cada caso, multiplicado pela probabilidade daquele caso ocorrer. $$ F[x] = \sum_i p[i]C[i] $$ Antes de começar, vejamos os casos extremos. Quando o vetor é vazio, não precisa de nenhuma comparação, então ##F[0]=0##. Se o vetor tiver tamanho um, então só precisa de uma comparação, ##F[1]=1##. É sempre bom ter uns casos pequenos para conferir no final.

Vamos ver agora o caso de um vetor de tamanho ##n##. Eu não sei qual número vai ser escolhido para cortar o vetor em dois, é o gerador de números aleatórios que escolhe. Então eu vou tirar a média sobre todas as escolhas possíveis. Como os ##n## valores são equiprováveis, então a probabilidade individual é ##1/n##: $$F[n] = \sum_{0\le k\lt n}\frac{1}{n}F[k,n] = \frac{1}{n}\sum_{0\le k\lt n}F[k,n] $$ Na fórmula acima, ##F[n]## é o número médio de comparações para um vetor de tamanho ##n##, e ##F[k,n]## é o número médio de comparações para o vetor de tamanho ##n##, no caso em que o corte foi feito no elemento ##k##.

Vamos calcular ##F[k,n]## agora. Esse valor depende de onde está o número que estamos procurando. Eu não sei onde está o número; como estamos calculando o caso médio, precisamos tirar a média sobre todos as posições possíveis. $$F[k,n]=\sum_{0\le i\lt n}\frac{1}{n}F[i,k,n]=\frac{1}{n}\sum_{0\le i\lt n}F[i,k,n]$$ Agora ##F[i,k,n]## significa o valor médio para um vetor de tamanho ##n##, cortado na posição ##k##, e com número final na posição ##i##.

Nós temos três casos agora. Se ##i=k##, então você só precisa comparar uma vez, a comparação vai ser verdadeira e o algoritmo termina. Se ##i\lt k##, então você compara uma vez e faz recursão para o lado esquerdo. Se ##i\gt k##, então você compara uma vez e faz recursão para o lado direito. Resumindo: $$F[i,k,n]=\begin{cases} 1+F[k] &(i \lt k) \\ 1 &(i = k) \\ 1+F[n-k-1] &(i \gt k) \\ \end{cases}$$ Podemos sintetizar tudo que vimos em uma única recorrência. $$\begin{align*} F[0] &= 0 \\ F[n] &= \frac{1}{n} \sum_{0\le k \lt n}\left( \frac{1}{n}\sum_{0\le i\lt k}\left(1+ F[k]\right) + \frac{1}{n} + \frac{1}{n}\sum_{k\lt i\lt n}\left(1+ F[n-k-1]\right) \right) \end{align*}$$ Agora é só resolver!

Antes de mais nada, vamos coletar quem é constante em relação ao índice da somatória. Todos os ##1/n## coletam, e as somatórias mais internas viram constantes. $$\begin{align*} F[n] &= \frac{1}{n^2} \sum_{0\le k \lt n}\left( \sum_{0\le i\lt k}\left(1+ F[k]\right) + 1 + \sum_{k\lt i\lt n}\left(1+ F[n-k-1]\right) \right) \\ &= \frac{1}{n^2} \sum_{0\le k \lt n} k \left(1+ F[k]\right) + 1 + \left(n-k-1\right)\left(1+ F[n-k-1]\right) \\ &= \frac{1}{n^2} \sum_{0\le k \lt n} k + k F[k] + 1 + n-k-1 +\left(n-k-1\right)F[n-k-1] \\ &= \frac{1}{n^2} \sum_{0\le k \lt n} n + k F[k] +\left(n-k-1\right)F[n-k-1] \\ &= 1+\frac{1}{n^2} \left(\sum_{0\le k \lt n} k F[k] +\sum_{0\le k \lt n} \left(n-k-1\right)F[n-k-1] \right) \\ \end{align*}$$ Vamos focar naquela segunda somatória. Primeiro fazemos ##q=n-k-1##: $$ \sum_{0\le k \lt n} \left(n-k-1\right)F[n-k-1] =\sum_{0\le n-q-1 \lt n} qF[q] $$ Agora é só manipular os índices: $$\begin{align*} 0\le n-&q-1 \lt n \\ 1-n\le -&q \lt 1 \\ -1\lt &q \le n-1 \\ 0\le &q \lt n \\ \end{align*}$$ Olha só, a segunda somatória é igual à primeira! $$ \sum_{0\le n-q-1 \lt n} qF[q] =\sum_{0\le q \lt n} qF[q] =\sum_{0\le k \lt n} kF[k] $$ Substituindo, chegamos em uma fórmula curta para ##F[n]##: $$\begin{align*} F[n] &= 1+\frac{1}{n^2} \left(\sum_{0\le k \lt n} k F[k] +\sum_{0\le k \lt n} \left(n-k-1\right)F[n-k-1] \right) \\ &= 1+\frac{1}{n^2} \left(\sum_{0\le k \lt n} k F[k] +\sum_{0\le k \lt n} kF[k] \right) \\ &= 1+\frac{2}{n^2} \sum_{0\le k \lt n} k F[k] \\ \end{align*}$$ A somatória ainda está atrapalhando, o ideal seria sumir com ela. Uma maneira de fazer isso é isolando: $$\begin{align*} F[n] &= 1+\frac{2}{n^2} \sum_{0\le k \lt n} k F[k] \\ \sum_{0\le k \lt n} k F[k] &= \frac{n^2}{2}\left(F[n]-1\right) \end{align*}$$ Essa última fórmula vale para todo ##n##, em especial vale também para ##n-1##: $$\begin{align*} \sum_{0\le k \lt n} k F[k] &= \frac{n^2}{2}\left(F[n]-1\right) \\ \sum_{0\le k \lt n-1} k F[k] &= \frac{(n-1)^2}{2}\left(F[n-1]-1\right) \\ \end{align*}$$ Ahá! Agora é só subtrair uma da outra que a somatória desaparece! $$\begin{align*} \sum_{0\le k \lt n} k F[k] - \sum_{0\le k \lt n-1} k F[k] &= \frac{n^2}{2}\left(F[n]-1\right) - \frac{(n-1)^2}{2}\left(F[n-1]-1\right)\\ (n-1) F[n-1] &= \frac{n^2}{2}\left(F[n]-1\right) - \frac{(n-1)^2}{2}\left(F[n-1]-1\right)\\ (n-1) F[n-1] &= \frac{n^2}{2}F[n]-\frac{n^2}{2} - \frac{(n-1)^2}{2}F[n-1]+\frac{(n-1)^2}{2}\\ \frac{n^2}{2}F[n] &= \left(\frac{(n-1)^2}{2}+(n-1)\right)F[n-1]+\left(\frac{n^2}{2} -\frac{(n-1)^2}{2}\right)\\ n^2F[n] &= \left(n^2-1\right)F[n-1]+(2n-1)\\ \end{align*}$$ Chegamos finalmente em uma recorrência sem somatórias. Melhor ainda, essa recorrência está no ponto certo para usar a técnica do summation factor! Como esse é um truque muito útil de se conhecer, eu vou fazer em câmera lenta. Você pode usar um summation factor sempre que sua recorrência for da forma abaixo: $$a_n F[n] = b_n F[n-1]+c_n$$ Esse é o nosso caso, se usarmos a correspondência abaixo: $$\begin{align*} a_n &= n^2 \\ b_n &= n^2-1 \\ c_n &= 2n-1 \end{align*}$$ O segredo do summation faction é tentar achar uma sequência mágica ##s_n## que tenha a seguinte propriedade: $$s_n b_n=s_{n-1}a_{n-1}$$ E por que isso é bom? Olha só o que acontece quando você multiplica a recorrência original por ##s_n##: $$\begin{align*} a_n F[n] &= b_n F[n-1]+c_n \\ s_n a_n F[n] &= s_n b_n F[n-1]+s_n c_n \\ s_n a_n F[n] &= s_{n-1}a_{n-1}F[n-1]+s_n c_n \\ \end{align*}$$ Agora, se você substituir ##T[n]=s_n a_n F[n]##, temos: $$\begin{align*} s_n a_n F[n] &= s_{n-1}a_{n-1}F[n-1]+s_n c_n \\ T[n] &= T[n-1]+s_n c_n \end{align*}$$ Opa, essa é fácil! Dá pra resolver de cabeça, é praticamente a definição de somatória! $$\begin{align*} T[n] &= T[n-1]+s_n c_n\\ T[n] &= T[0] + \sum_{1\le k \le n} s_k c_k \\ s_n a_n F[n] &= s_0 a_0 F[0] + \sum_{1\le k \le n} s_k c_k \\ F[n] &= \frac{1}{s_n a_n} \left(s_0 a_0 F[0] + \sum_{1\le k \le n} s_k c_k \right) \\ \end{align*}$$ Pronto, agora não é mais uma recorrência, é uma fórmula simples. A dificuldade do método é encontrar o tal do ##s_n##. Para achá-lo você pode usar a intuição, chutar, fazer um sacrifício aos deuses; ou então você pode usar o método do porradão. Nesse método você começa com ##s_n=s_{n-1}\left(a_{n-1}/b_n\right)## e abre recursivamente, chegando no monstrinho abaixo: $$s_n = \frac{a_{n-1}a_{n-2}\ldots a_2 a_1}{b_n b_{n-1}\ldots b_3 b_2}$$ Essa fórmula parece grande, e é grande mesmo. Mas na maioria das vezes tem um monte de cancelamentos, o que deixa o resultado final pequeno. Vamos ver no nosso caso como fica: $$\begin{align*} s_n &= \frac{a_{n-1}a_{n-2}\ldots a_2 a_1}{b_n b_{n-1}\ldots b_3 b_2} \\ &= \frac{(n-1)^2 (n-2)^2 \ldots 2^2 1^2} {\left(n^2-1\right) \left(\left(n-1\right)^2-1\right) \ldots (3^2-1) (2^2-1)} \\ &= \frac{(n-1)^2 (n-2)^2 \ldots 2^2 1^2} {\left((n+1)(n-1)\right) \left((n-1+1)(n-1-1)\right) \ldots (3+1)(3-1) (2+1)(2-1)} \\ &= \frac{(n-1) (n-2)(n-3) \ldots3\cdot 2\cdot 1} {(n+1) (n) (n-1)\ldots (4+1) (3+1) (2+1)} \\ &= \frac{2} {(n+1) (n) } \\ \end{align*} $$ Cancelou mesmo! Bastou lembrar que ##(a^2-b^2)=(a+b)(a-b)## que o resto saiu. Agora é só lembrar que ##F[0]=0## e substituir na fórmula: $$\begin{align*} F[n] &= \frac{1}{s_n a_n} \left(s_0 a_0 F[0] + \sum_{1\le k \le n} s_k c_k \right) \\ F[n] &= \frac{n(n+1)}{2n^2} \sum_{1\le k \le n} \frac{2(2k-1)}{k(k+1)}\\ \end{align*}$$ Essa é a solução da recorrência. Dá para achar uma forma fechada? Nesse caso dá, desde que você tope uma forma fechada em função dos números harmônicos ##H[n]##, que são definidos como: $$H[n]=\sum_{1\le k \le n}\frac{1}{k}$$ Você começa abrindo a fração no somatório: $$\begin{align*} \frac{2(2k-1)}{k(k+1)} &= \frac{A}{k} + \frac{B}{k+1}\\ \frac{4k-2}{k(k+1)} &= \frac{A(k+1)+B k}{k(k+1)}\\ 4k-2 &= Ak+A +B k \\ 4k-2 &= (A+B)k+A \\ A &= -2 \\ B &= 6 \end{align*}$$ Agora você divide o somatório em dois: $$\begin{align*} F[n] &= \frac{n(n+1)}{2n^2} \sum_{1\le k \le n} \frac{2(2k-1)}{k(k+1)}\\ F[n] &= \frac{n(n+1)}{2n^2} \sum_{1\le k \le n} \frac{6}{k+1} - \frac{2}{k}\\ F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - \sum_{1\le k \le n} \frac{2}{k} \right) \\ F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - 2\sum_{1\le k \le n} \frac{1}{k} \right) \\ F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - 2H[n] \right)\\ \end{align*}$$ O primeiro somatório precisa de um pouco de massagem: $$\begin{align*} \sum_{1\le k \le n} \frac{6}{k+1} &= \sum_{1\le q-1 \le n} \frac{6}{q} \\ &= \sum_{2\le q \le n+1} \frac{6}{q} \\ &= -6+\frac{6}{n+1}+6\sum_{1\le q \le n} \frac{1}{q} \\ &= -6+\frac{6}{n+1}+6H[n] \\ \end{align*}$$ Por fim: $$\begin{align*} F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - 2H[n] \right)\\ F[n] &= \frac{n(n+1)}{2n^2} \left( -6+\frac{6}{n+1}+6H[n] - 2H[n] \right)\\ F[n] &= -3 +\frac{2(n+1)}{n}H[n]\\ \end{align*}$$ Ufa, deu trabalho! E isso é ##O(\log n)## mesmo? Bem, podemos verificar usando um assintótico para o harmônico, por exemplo: $$ H[n] = \ln n + \gamma + O(1/n) $$ Agora você substitui: $$\begin{align*} F[n] &= -3 +\frac{2(n+1)}{n}H[n]\\ &= -3 +2H[n] +\frac{H[n]}{n} \\ &= -3 +2\ln n +2\gamma +O(1/n)+ 2\frac{\ln n}{n}+2\frac{\gamma}{n}+O(1/n^2) \\ &= 2\ln n + O\left(\frac{\ln n}{n}\right) \\ &= O(\log n) \end{align*}$$ Como queríamos demonstrar!

Ou seja, realmente a Busca por Bozo Binário é ##O(\log n)##. Em demonstrações grandes assim eu sempre gosto de fazer uma simulação por Monte Carlo para conferir as contas. A minha simulação está no github:

Busca por Bozo Binário, simulação por Monte Carlo

Simulando 50000 vezes em um vetor de tamanho 3000, o número médio de comparações foi de ##14.1769##. O valor teórico previsto pela fórmula é de ##14.1732##, nada mau!

Eu também fiz uma simulação da busca binária tradicional como comparação. Com os mesmos parâmetros, a busca tradicional usou ##10.6356## comparações no caso médio, ou seja, foi mais eficiente que o Bozo Binário.

Por que isso acontece? Analisando o assintótico fica claro. Para ##n## grande, o valor médio da Busca por Bozo é de ##2 \ln n## (logaritmo natural), enquanto que na busca binária é ##\log_2 n## (logaritmo de base 2). Então, em média, esperamos que Bozo Binário seja mais lento por um fator de:

$$\begin{align*} \frac{2\ln n}{\log_2 n} &= {2\frac{\log n}{\log e}} \div {\frac{\log n}{\log 2}} \\ &= 2\frac{\log n}{\log e} \times \frac{\log 2}{\log n} \\ &= 2\frac{\log 2}{\log e} \\ &= 2\ln 2 \\ &\simeq 1.386 \end{align*}$$ Ou seja, mais ou menos uns 38% mais lento, o que bate com a simulação.

Depois da análise do algoritmo, a dúvida que resta é: e o candidato? Bem, se fosse eu a entrevistar, precisaria coletar mais dados para saber se ele realmente sabia do tempo médio, ou se estava perdido e implementou a primeira coisa que veio na cabeça.

Um possível follow-up é "Certo, esse método é mais complicado que a busca binária tradicional, mas funciona. Você pode me falar alguma situação onde ele é preferível sobre o tradicional?". Curiosamente, existe sim um caso onde a Busca por Bozo Binário é melhor que a busca binária tradicional! Mas esse caso fica para a próxima :)

A Intuição do Knuth

Às vezes eu me pergunto se as pessoas da minha área têm noção de quão sortudos nós somos. Os físicos adorariam viajar no tempo para conversar com o Newton, os matemáticos adorariam conversar com o Euclides, os biólogos adorariam conversar com o Darwin. Mas nós podemos conversar com o Knuth!

Nós temos a sorte de viver no mesmo período de tempo que o criador da análise de algoritmos, que é uma das bases da Ciência da Computação. Se você gosta do assunto, vale a pena juntar uns trocos e viajar até a Califórnia para assistir a uma das palestras dele (dica: todo fim de ano, inspirado nas árvores de Natal, ele faz uma palestra de estrutura de dados, falando sobre árvores; elas também estão online se você não tiver como ver ao vivo).

Eu fiz a peregrinação em 2011, quando consegui assistir a uma das palestras dele. Aproveitei para ir todo contente pegar minha recompensa por ter achado um erro no Art of Computer Programming, mas ele, marotamente, me disse que aquilo que eu achei não era um erro, era uma pegadinha, e eu caí! (Mas eu não vou falar qual a pegadinha, vá na página 492 do TAOCP volume 4A, primeira edição, e confira você mesmo :)

Eu e Knuth, o trollzinho

Nesse dia perguntaram que opinião ele tinha sobre o problema mais difícil da nossa geração, P=NP. A intuição dele é que provalmente é verdade, mas ele acredita que se acharmos a demonstração, ela vai ser não-construtiva. O que isso significa? O que é uma demonstração não-construtiva?

Demonstrações construtivas e não-construtivas

Em análise de algoritmos, as demonstrações construtivas são as mais comuns. Por exemplo, digamos que eu quero provar que é possível calcular x elevado a y em tempo O(y). Isso é fácil, basta construir um algoritmo assim:
E se eu quiser provar que esse mesmo problema pode ser resolvido em tempo O(log y)? Novamente, tudo que eu preciso fazer é exibir um algoritmo que implemente isso:

(Nesse caso eu também precisaria provar que esse algoritmo é de fato O(log y), já não é óbvio por inspeção). Nos dois casos temos exemplos de demonstrações construtivas: se eu quero provar uma propriedade P, basta exibir um algoritmo que tenha essa propriedade P.

As demonstrações não-construtivas são diferentes. Nelas, eu posso provar a propriedade P sem mostrar o algoritmo, através de alguma propriedade matemática do modelo.

Por exemplo, imagine que eu tenho uma lista ordenada de números. Se eu fizer uma busca binária, posso achar a posição de um número dado com O(log n) comparações. Mas é possível criar um algoritmo mais rápido que isso? Eu digo que não é possível, e para isso vou fazer uma prova não-construtiva de que esse é o mínimo que um algoritmo de busca precisa para funcionar.

A Teoria de Shannon aplicada à busca binária

Para isso eu vou usar a teoria da informação de Shannon. Essa teoria é surpreendentemente intuitiva, e se baseia no conceito de surpresa. Se eu te falar que o céu ficou escuro às 19h, você não vai achar nada de mais, nessa hora o Sol está se pondo, então é natural que o céu fique escuro. Mas e se eu falar que o céu ficou escuro às 10 da manhã? Foi uma tempestade? Um eclipse? A nave do Independence Day?

Intuitivamente, quanto mais surpresos nós ficamos com uma sentença, mais informação ela tem. O Shannon definiu então a quantidade de informação como sendo uma função monotônica da probabilidade do evento acontecer:

Se o evento é raro, tem bastante informação; se o evento é comum, tem pouca informação. A base do logaritmo fornece a unidade de medida, se a base for 2, então a informação é medida em bits.

E quanta informação nós ganhamos com uma comparação? Se a chance de dar verdadeiro ou falso for a mesma, então a chance é p(m)=1/2, logo a informação é I(m)=1. Você ganha exatamente um bit de informação com uma comparação.

Qual o resultado do nosso algoritmo de busca? O resultado é um índice, se nós temos n elementos no vetor, então a resposta é um índice que varia de 0 a n-1. Logo, a probabilidade de você escolher o índice certo ao acaso é p(m)=1/n, já que a escolha é uniforme.

Quanta informação tem essa escolha, então? Fazendo a conta:

Se você precisa de log n bits para descrever a resposta, e você ganha só 1 bit por comparação, então não tem como um algoritmo rodar em menos que O(log n): a informação tem que vir de algum lugar! Com isso, nós mostramos que qualquer algoritmo precisa rodar no mínimo em tempo O(log n), e sem precisar mostrar o algoritmo em si. Essa é uma demonstração não-construtiva.

Pressinto a pergunta: "mas RicBit, e a busca com hash table, ela não é O(1)?". Sim, ela é! Mas ela não usa comparações, e a nossa análise foi exclusivamente para métodos baseados em comparações. Com um acesso a uma hash você pode ganhar mais que 1 bit de informação por operação.

O limite da ordenação

Um outro exemplo é achar o limite dos algoritmos de ordenação. Suponha que eu tenho um vetor com elementos bagunçados e quero ordená-los usando comparações. Eu sei que cada comparação ganha 1 bit de informação, então só preciso saber quanta informação tem na saída.

Qual o resultado do algoritmo? Um vetor ordenado. Mas os valores do vetor em si são irrelevantes, o que importa mesmo é saber a ordem relativa entre eles. Essa ordem relativa pode ser expressa como uma permutação dos itens originais.

Quantas permutações existem? Se o vetor tem tamanho n, então existem n! permutações, logo a probabilidade é 1/n!. Fazendo as contas:

Primeiro você usa a aproximação de Stirling, depois joga fora todos os termos assintoticamentes menores que o dominante. O resultado é que nós provamos que nenhuma ordenação pode ser melhor que O(n log n), sem precisar mostrar nenhum algoritmo!

Novamente, esse resultado só vale para ordenações baseadas em comparações. Sem usar comparações, você tem métodos como radix sort e ábaco sort que são melhores que O(n log n).

A análise por quantidade de informação

Esse método de análise da quantidade de informação pode ser utilizado em qualquer algoritmo, desde que você note um detalhe muito importante: o método acha um limite inferior para a complexidade, mas não prova que esse algoritmo existe! Tudo que conseguimos provar como ele é que, se o algoritmo existir, então ele não pode ser melhor que o limite achado.

O Jogo do Pi

No último post eu falava de variações em jogos, aí lembrei de uma variação divertida em um jogo bem conhecido. Quem assistia Topa Tudo Por Dinheiro certamente deve se lembrar do Jogo do Pim. O objetivo é enumerar os naturais pelo maior tempo possível, trocando os múltiplos de quatro pela palavra "pim":

1,2,3,pim,  5,6,7,pim, 9,10,11,pim...

Parece simples, mas na prática muita gente se confundia e errava antes mesmo de chegar ao quarenta!

Uma variação mais difícil desse jogo foi inventada pelo Juca uns anos atrás: o Jogo do Pi. Dessa vez você ainda precisa enumerar os naturais, mas precisa falar "π" toda vez que passar um múltiplo inteiro de pi. Os primeiros múltiplos inteiros são aproximadamente 3.14, 6.28 e 9.42, então a sequência começa assim:

1,2,3,π,  4,5,6,π, 7,8,9,π...

Parece o mesmo jogo né? Você conta três números, fala π, conta mais três, fala π, e assim por diante. Mas se você seguir essa estratégia, vai perder! Olha o que acontece se você continuar:

..., 16,17,18,π, 19,20,21,π, 22,23,24,25,π,...

A estratégia de contar de três em três falha! Entre 7π e 8π tem quatro naturais ao invés de só três.

Bem, será que não dá pra melhorar a estratégia? Podemos contar quantos números tem entre os múltiplos de pi. Com sorte, o padrão é 3‑3‑3‑3‑3‑3‑3‑4 e depois repete. Nesse caso, o padrão do Jogo do Pi é oito vezes maior que o padrão do Jogo do Pim. Infelizmente, se você fizer as contas, esse padrão também é quebrado após 112π.

Será que devemos procurar um padrão de padrões então? Nessa altura não já compensa fazer as contas na mão, melhor deixar o python achar o período pra gente:

Script em python para tentar achar um período

Más notícias: o script acha padrões maiores e maiores, sem parar. Não tem como concluir se existe ou não um padrão só analisando a saída do programinha.

O jeito de resolver essa dúvida, então, é usando matemática! Se você tem medo de teoria dos números, pule o quadro azul:

Antes de mais nada, vamos colocar nosso problema na notação correta. Nós queremos descobrir se a quantidade de números entre múltiplos inteiros consecutivos de pi formam uma sequência periódica. Vamos primeiro escrever a sequência explicitamente:

Isso gera a sequência que queremos. Agora vamos supor que essa sequência tem um período p. Das duas uma: ou a gente acha o valor de p, ou batemos em uma contradição pelo caminho.

Se nós somarmos todos os primeiros n termos dessa sequência, temos uma soma telescópica e um resultado curto:

Isso vale para todo n, inclusive no caso onde n é um múltiplo inteiro do período p:

Mas, nesse caso, podemos fazer a somatória de outro jeito. Ao invés de somar direto de 0 a kp, podemos somar k vezes o período p. Aí, como a sequência é periódica, um dos termos some:

O termo mais interno da somatória não varia mais com i, então essa somatória está somando k vezes a mesma coisa:

A conclusão é que, quando p é o período da sequência, podemos jogar um k inteiro de dentro para fora do piso. Podemos agora abrir a definição de piso para kpπ:

Dividindo todo mundo por kp, temos:

Agora, para todo real ε > 0, podemos escolher um k grande o suficiente tal que ε seja maior que 1/kp. Logo:

Pelo teorema do sanduíche, concluímos que:

Note que o numerador dessa função é um piso, logo é um inteiro. Sabemos também que p é um inteiro, logo a fração é um racional. Mas pi é irracional, portanto chegamos a uma contradição, e a sequência não tem período.

Ou seja, concluímos matematicamente que o Jogo do Pi não tem período. Outra maneira de dizer a mesma coisa é que o período do Jogo do Pi é infinito, e portanto o Jogo do Pi é infinitamente mais díficil que o Jogo do Pim!

Agradecimentos ao povo do Math Stack Exchange pela ajuda na demonstração

A Fórmula de Bhaskara

Um tempo atrás o Karlisson teve uma idéia muito boa: criar uma coletânea com os 1001 algoritmos que você deve implementar antes de morrer. O objetivo é fazer uma implementação de referência, em python, para todos os algoritmos clássicos, desde os simples como o Bubblesort até os complexos como o Edmonds-Karp.

Se você também quiser participar, é só fazer um fork do repositório no github. É um ótimo exercício para quem é iniciante, e para os veteranos é uma boa oportunidade de praticar o hábito de fazer code reviews. Já que no post anterior eu falava de parábolas, vamos aproveitar pra revisar um dos códigos do projeto: uma função que acha as raízes reais da equação de segundo grau:

def bhaskara(a, b, c):
    delta = b ** 2 - 4 * a * c
    if delta < 0:
        return None
    else:
        raizes = []
        m1 = math.sqrt(delta)
        r1 =(-b + m1) / (2 * a)
        raizes.append(r1)
        r2 =(-b - m1) / (2 * a)
        raizes.append(r2)
        return raizes

A função é uma implementação direta da fórmula:

Conferindo com a fórmula, o código parece correto. Porém, olhando com cuidado, ele tem três problemas!

Primeiro problema: API inconsistente

O problema mais imediato do código é que ele faz uma divisão por zero quando a é igual a zero. Mas isso é um bug de verdade? A intenção do código era resolver equações do segundo grau, mas, se a é nulo, então a equação não é de segundo grau, é de primeiro!

Eu interpreto isso como um problema na API. Se a idéia do código é resolver somente equações de grau 2, então, antes de sair fazendo conta, ele deveria verificar se as entradas são válidas (ou seja, se a é não-nulo). Por outro lado, se a intenção é resolver equações de grau 2 ou inferior, então ele poderia retornar [-c/b] quando a for zero. Nos dois casos, faltou sinalizar exatamente qual o propósito da API.

A API ainda tem uma segunda falha, que é o valor de retorno. O que exatamente ele vai retornar? Eu poderia, por exemplo, dizer que a função sempre retorna uma lista com as soluções existentes. Nesse caso, quando o delta é negativo, o correto seria retornar uma lista vazia, ao invés de retornar None.

Além disso, o que devemos retornar quando o delta é exatamente zero? As duas possibilidades são retornar uma lista com uma única raiz, [-b/2a], ou então retornar uma lista com duas raízes iguais, ou seja, [-b/2a,-b/2a]. Qual resposta faz mais sentido?

A escolha correta vem do Teorema Fundamental da Álgebra, que diz que o número de raízes complexas de uma equação polinomial é igual ao grau da equação. O motivo de considerarmos que algumas equações tem múltiplas raízes iguais é pra fazer esse teorema ser válido em todas as situações.

Mas o teorema só vale para raízes complexas! Como estamos trabalhando com raízes reais, então eu acho que o mais correto é retornar uma lista com um único elemento mesmo.

Você também poderia argumentar que o nome da função não é apropriado, porque está descrevendo a implementação (fórmula de Bhaskara) ao invés de descrever a interface (solução da equação quadrática). Mas aí nós chegamos no segundo problema...

Segundo problema: Bhaskara quem?

Certamente você ouviu na escola que a solução da equação quadrática chama-se fórmula de Bhaskara. Mas, curiosamente, é só no Brasil que a fórmula tem esse nome!

Pode conferir: a wikipedia em inglês nem cita o Bhaskara. Talvez na wikipedia em hindi? Nada de Bhaskara lá também, na Índia eles chamam a equação de fórmula de Sridhar Acharya. Nos outros países a fórmula nem tem nome, é simplesmente "equação quadrática" ou "fórmula a-b-c".

Na verdade, os babilônicos já sabiam resolver alguns tipos de equações quadráticas há pelo menos 4000 anos. Os gregos conheciam soluções geométricas. E a fórmula em si, é do Bhaskara mesmo?

Bem, uma maneira de resolver a dúvida é procurando exatamente o que ele escreveu. O texto mais importante do Bhaskara foi o livro que dedicou à filha, Lilavati. O original era em sânscrito, mas tem uma tradução em inglês online para quem quiser ler. Eu procurei pela solução da quadrática, e o mais próximo que achei foi a proposição 63:

Uma quantidade, incrementada ou decrementada de sua raiz quadrada multiplicada por algum número, é dada. Então adicione o quadrado de metade do multiplicador da raiz ao número dado, e extraia a raiz quadrada da soma. Adicione metade do multiplicador, se a diferença foi dada; ou subtraia, se a soma for dada. O quadrado do resultado é a quantidade procurada.

Em notação moderna, o que o Bhaskara escreveu foi:

Eu não acho que seja correto creditar ao Bhaskara a solução da equação quadrática. A fórmula acima é bem parecida com a que usamos, mas note o detalhe: ela só fornece uma raiz! Para considerar a solução como correta, eu acho que ele deveria ter explicitado que algumas equações vão ter duas raízes.

Pois bem, podemos então mudar o nome da função, de bhaskara() para algo mais apropriado, como quadratic_roots(). Mas ainda resta o último problema...

Terceiro problema: Instabilidade numérica

Vamos testar o código original com a seguinte equação: x^2-2000000000x+4=0. As raízes são aproxidamadamente 2e9 e 2e-9. Porém, olhe só o resultado:

>>> bhaskara.bhaskara(1, -2e9, 4)
[2000000000.0, 0.0]

Como assim zero? Cadê a solução 2e-9? Essa raiz foi vítima de um problema que nem todo mundo presta atenção: muito cuidado com a perda de precisão em ponto flutuante. Em especial, quando você subtrai dois números de magnitude similar, a resposta pode não ter precisão suficiente para dar o resultado correto. Para a equação dada, o código original tem os valores b e m1 muito próximos, e a subtração perde a menor raiz.

A solução é evitar a subtração. Você inspeciona o sinal de b, e escolhe m1 com o mesmo sinal, de modo que os módulos são sempre somados, achando assim a primeira raiz x1. Tendo x1, você pode achar x2 usando a fórmula de Viète para o produto de raízes:

Dessa maneira a subtração problemática some, e chegamos na versão final do programa corrigido:

def quadratic_roots(a, b, c):
  """Find all real roots of a quadratic equation.

  Args:
    a, b, c: Coefficients of the quadratic, a*x*x+b*x+c=0
  Returns:
    A list with all real roots.
  Raises:
    ValueError if a==0 (the equation is not quadratic).
  """
  if not a:
    raise ValueError
  delta = b * b - 4 * a * c
  if delta < 0:
    return []
  if not delta:
    return [-b / (2 * a)]
  if b < 0:
    x1 = (-b + delta ** 0.5) / (2 * a)
  else:
    x1 = (-b - delta ** 0.5) / (2 * a)
  x2 = c / (a * x1)
  return [x1, x2]

Conferindo:

>>> roots.quadratic_roots(1, -2e9, 4)
[2000000000.0, 2.0000000000000001e-09]

Nice!

Você pode se perguntar se esse exemplo que eu dei não é artificial demais, e se coisas assim acontecem na prática. Mas acontecem sim! Um exemplo comum onde esse erro acontece é na implementação da busca binária.

Por exemplo, suponha que você quer resolver a equação ln(x)+sqrt(x)=5. A função ln() é monotônica, o sqrt() também, então dá pra achar a raiz por busca binária:

def search(epsilon):
  inf, sup = 0.0, 1e10
  while sup - inf > epsilon:
    med = (sup + inf) / 2
    x = math.log(med) + math.sqrt(med)
    if x > 5.0:
      sup = med
    else:
      inf = med
  return inf

Vamos testar essa rotina para várias precisões:

>>> binary.search(1e-4)
8.3093709690729156
>>> binary.search(1e-10)
8.3094326941933723
>>> binary.search(1e-15)

Oops! Quando epsilon vale 1e-15, a rotina trava! De novo, o problema é perda de precisão. Nesse caso, as variáveis sup e inf tem valores muito próximos, então a média med não sai do lugar, ficando presa no valor do inf para sempre.

A solução, novamente, é sumir com a subtração. Ao invés de controlar os extremos do intervalo, você usa uma variável para o valor inicial e outra para o tamanho do intervalo:

def smart_search(epsilon):
  inf, delta = 0.0, 1e10
  while delta > epsilon:
    delta /= 2
    med = inf + delta
    x = math.log(med) + math.sqrt(med)
    if x < 5.0:
      inf += delta
  return inf

Conferindo:

>>> binary.smart_search(1e-4)
8.3093709690729156
>>> binary.smart_search(1e-10)
8.309432694193374
>>> binary.smart_search(1e-15)
8.3094326942315693
>>> binary.smart_search(1e-30)
8.3094326942315693

Agora sim!

No fim das contas, é importante você ter a matemática sempre afiada, mas também é igualmente importante conhecer os limites impostos pela computação do mundo real.

(Agradecimentos ao Karlisson, ao Henrique, ao Wladimir, ao Chester e ao povo da lista eng-misc pela ajuda na pesquisa!)

Sorteios e Aniversários

No fim de semana passado eu participei do IV ENSL: o Encontro Nordestino de Software Livre, que nessa última edição foi realizado em Natal, no Rio Grande do Norte. Eu achei bastante apropriado, dado que Natal sempre foi um pólo de inovação tecnólogica! Citando alguns exemplos:

• Natal é a sede do Instituto de Neurociências do Miguel Nicolelis, pioneiro das interfaces cérebro-máquina e um dos maiores cientistas vivos.
• Alex Kipman, diretor de incubação da Microsoft, se inspirou na cidade para criar o Xbox Kinect (que, vocês devem lembrar, tinha o nome-código de Project Natal).
• Na década de 80, era onde estudava Tadeu Curinga da Silva, pioneiro dos videogames nacionais e autor do jogo Em Busca dos Tesouros, o melhor jogo da plataforma ZX-81.
• Além disso, Natal também é onde mora o Karlisson Bezerra, autor das tirinhas do Nerdson!

Embora o turismo nerd tenha sido excelente, eu fiquei um pouco frustrado de não ter aproveitado tanto o turismo tradicional. Eu consegui ir à praia e provar a deliciosa manteiga de garrafa, mas não deu tempo de fazer o típico passeio de buggy pelas dunas. Por outro lado, posso não ter feito o passeio de buggy, mas presenciei um outro tipo de bug!

Tradicionalmente, todo encontro de software livre termina com um sorteio de brindes. Como é um evento de programadores, naturalmente não se usou um saco plástico cheio de papéis com os nomes de cada um. Ao invés disso, o povo programou uma solução na hora mesmo, usando só três linhas de python:

import random
x = open("participantes.txt").readlines()
print x[random.randint(0, len(x))]

Como isso foi feito no modo interativo, então basta repetir a última linha para realizar um novo sorteio. Mas enquanto eles digitavam isso no telão, eu pensei: putz, esse código tem dois bugs. Você consegue achá-los?

Primeiro Bug

O primeiro bug, na real, não é do programa e sim da API do Python. Quando você projeta uma API, precisa ter muito cuidado para que as interfaces sejam consistentes, e o randint() é um dos raros casos onde o Python pisa no tomate.

O problema é que quase todas as definições de intervalo do Python usam intervalos semi-abertos, mas o randint usa um intervalo fechado (se você não lembra o que é um intervalo aberto ou fechado, pergunte ao professor Joe). Na notação de list slices, ou mesmo no range(), o último número não faz parte da lista; mas no randint() o último número conta.

Por isso, se você tentar sortear com randint(0, len(x)), ele pode sortear um número além do fim da lista, e você vai ter um erro do tipo list index out of range. Murphy sendo Murphy, é claro que isso aconteceu ao vivo lá no telão :)

Depois que uma API já está publicada e sendo usada, é muito difícil consertar um bug assim. Num mundo ideal você consertaria a interface, mas isso poderia quebrar todos os códigos legacy que a usam. Por isso, o melhor que dá pra fazer é depreciar a função original e propor uma interface nova, que no nosso caso é a randrange():

print x[random.randrange(0, len(x))]

Nesse problema em específico tem uma outra solução até mais elegante, que é usar o choice():

print random.choice(x)

Mas nenhuma dessas soluções conserta o outro bug do programa...

Segundo bug

Esse é mais sutil. No começo o programa estava funcionando bem, mas em um certo ponto, eis que o programa sorteia uma pessoa que já tinha sido sorteada antes! Eram mais ou menos 700 participantes, então a chance disso acontecer deveria ser bem pequena né? Mas algum tempo depois dessa pessoa, teve mais uma que foi sorteada duas vezes, e depois dessa ainda mais uma. Como pode um evento tão raro acontecer tantas vezes?

Marmelada não era, porque todo mundo viu o código, então o erro deve estar no gerador de números aleatórios do Python, certo? Surpreendentemente, a resposta é não! Esse é um caso clássico onde a nossa intuição falha. A análise ingênua é que a chance de um número ser sorteado duas vezes é 1/700 vezes 1/700, ou seja, 0.0002%. Mas esta é a análise correta do problema errado.

Se eu tivesse feito só dois sorteios, entre 700 participantes, e perguntasse qual era a chance do Ricbit ganhar nos dois, aí sim a chance seria 0.0002%. Mas o nosso problema é diferente. Nós temos um número p=700 de participantes, e vamos fazer s=70 sorteios (de olho eu acho que foram uns 70 sorteios: tinha brinde pra caramba no evento, e alguns sorteados eram descartados porque não estavam na sala).

O enunciado correto então é: qual o valor esperado do número de pessoas que é sorteada duas vezes? Esse valor é bem mais alto do que a intuição imagina!

Vamos calcular passo a passo quanto dá isso. Primeiro, vamos calcular qual é a probabilidade do Ricbit ser sorteado no primeiro evento. Isso é fácil né, a chance é 1/p. A chance oposta, ou seja, do Ricbit não ser sorteado no primeiro evento, é de (1-1/p). Continuando, a chance do Ricbit não ser sorteado em nenhum evento, é a abaixo:



Qual é a chance do Ricbit ser sorteado no primeiro evento e não ser sorteado em nenhum outro? Pelo mesmo raciocínio acima, ela é:



E qual a chance do Ricbit ser sorteado uma única vez durante toda a série de sorteios? Ora, é a chance de ser sorteado só no primeiro sorteio, mais a chance de ser sorteado só no segundo, e assim por diante. Como essas chances são todas iguais, então a chance do Ricbit ser sorteado uma única vez é:



Uma outra maneira de pensar a fórmula acima é que ela representa a chance de ser sorteado uma única vez numa posição qualquer, vezes o número de posições possíveis.

Agora chegamos onde queríamos, qual a chance do Ricbit ser sorteado exatamente duas vezes na série toda? Analogamente, é a chance de ser sorteado duas vezes, multiplicada pelo número de combinações possíveis dos dois sorteios ao longo da série. Logo, ela vale:



Por fim, qual o valor esperado do número de repetições? Uma repetição é quando uma pessoa é sorteada duas vezes: então basta, para cada pessoa, somar a probabilidade de ser sorteado duas vezes. Daí, o valor esperado final é:



Substituindo os parâmetros, descobrimos que o valor esperado é 3.1, o que realmente é bem próximo das três repetições que eu presenciei (nos cálculos eu desprezei sorteios triplos, quádruplos, etc, porque esses são realmente raros). Eu fiz uma simulação de Monte Carlo para quem só acredita vendo:

Simulação de Monte Carlo em Python

Na verdade, esse resultado contra-intuitivo é bem conhecido na literatura, o povo chama de Paradoxo do Aniversário. Numa festa com 30 pessoas, qual a chance de haver duas pessoas com o mesmo aniversário? A intuição é que isso deve ser raro, mas a probabilidade, na verdade, é de 70%. A argumentação é exatamente a mesma anterior, e, de fato, se você usar a nossa fórmula com s=30 e p=365, vai descobrir que o valor esperado é 1.1. Ou seja, em média, festas com 30 pessoas costumam ter um par de pessoas com o mesmo aniversário.

Se você trabalha com computação, então isso tem uma aplicação prática. A chance de repetição em um sorteio é exatamente a mesma chance de uma colisão em uma hash table, assumindo que sua função de hash tem distribuição uniforme. Por isso, quando você for dimensionar o tamanho de uma hash, pode usar nossa fórmula para saber o número esperado de colisões (por exemplo, em 70 inserções numa hash de tamanho 700, você vai ter em média 3 colisões).

E como resolver o problema original do sorteio? Basta tirar da lista os números que já foram sorteados! Uma implementação curta em python é a abaixo:

import random
x = open("participantes.txt").readlines()
random.shuffle(x)

print x.pop()

Ao invés de sortear nomes de maneira independente, você simplesmente faz uma permutação aleatória da lista com o shuffle() e vai tirando os elementos de um em um usando o pop(). Esse programa é bem curto e não tem os problemas do original. Mas, no fim das contas, se o programa original não tivesse nenhum bug, o sorteio não teria sido engraçado como foi :)

Se você é veterano do FISL, considere no ano que vem ir também ao ENSL. O povo é divertido, a comida local é excelente, e nada supera ir a uma praia do Nordeste entre uma palestra e outra :)

O algoritmo mais lento do oeste

Por que as pessoas aprendem complexidade computacional? Para achar os algoritmos mais rápidos, é claro. Dada uma lista de algoritmos equivalentes, você pode ordená-los pela complexidade, e o primeiro da lista é o mais eficiente.

Uma brincadeira divertida é fazer isso ao contrário: dada uma especificação, qual o algoritmo determinístico mais lento que a resolve? Por exemplo, qual o algoritmo mais lento para decompor um número em fatores primos?

O algoritmo mais rápido que resolve esse problema é o algoritmo de Shor, que roda em O((log n)³) e requer um computador quântico. Meu candidato a algoritmo mais lento é o abaixo, que na falta de nome melhor eu chamei de Fatoração Bozo:

import itertools

def factor(n):
 solutions = []
 for f in itertools.product(range(1,1+n),repeat=n):
   if reduce(lambda x,y: x*y, f) == n:
     solutions.append(filter(lambda x:x>1, list(f)))
 solutions.sort(key=len, reverse=True)
 return solutions[0]

Parece complexo, mas o funcionamento é bem simples. Vamos pensar no nosso objetivo: queremos achar uma lista de fatores primos, que multiplicados resultam no número desejado (vamos chamá-lo de N).

Nós sabemos que os primos dessa lista são menores ou iguais a N, e sabemos também que essa lista nunca vai ter mais que N elementos. A primeira coisa que fazemos no algoritmo é gerar todas as possíveis listas candidatas, que são o produto cartesiano do intervalo de 1 a N com N repetições. Por exemplo, para N=3, teríamos 27 combinações:

111 112 113 121 122 123 131 132 133
211 212 213 221 222 223 231 232 233
311 312 313 321 322 323 331 332 333

Em python o jeito mais fácil de gerar essa lista é usando itertools.product(), que faz o serviço sozinho. Mas, dessas listas, só interessam aquelas cujo produto seja igual a N. Por exemplo, se N=4, então as listas interessantes são:

1114 1122 1141 1212 1221
1411 2112 2121 2211 4111

Python não tem um operador para multiplicação de listas como tem o sum() para adições, mas pra isso que serve o reduce() né. As listas selecionadas até agora estão boas, mas elas têm um monte de elementos 1 que não contribuem em nada. Usando o filter() podemos jogar fora os uns. Para N=4:

4 22 4 22 22 4 22 22 22 4

Várias listas repetidas, mas tudo bem. O truque agora é notar que, de todas essas listas, a fatoração em primos é a lista de comprimento máximo. É fácil provar isso: se a lista máxima tivesse um elemento composto, então você poderia fatorá-lo, e trocar o número composto pela sua fatoração. Logo, você estaria gerando uma lista maior que a lista máxima, o que é uma contradição.

Pedindo para que o Python ordene as listas por tamanho, basta pegar a primeira lista da seqüência e temos então a nossa fatoração!

E qual a complexidade dessa solução? Ela é dominada pelo número de listas que geramos no primeiro passo, então a complexidade da Fatoração Bozo é maior que O(nⁿ). Se você não tem noção de como isso é lento, veja o tempo que esse código leva para fatorar alguns números:

N=6 - 0.2 segundos
N=7 - 2.1 segundos
N=8 - 45.8 segundos
N=9 - 20 minutos

A função cresce tão rapidamente que, para fatorar N=18, você levaria mais tempo que a idade do universo até agora! Como exercício para o leitor, tente calcular quanto tempo levaria para fatorar o RSA-704, que em decimal tem 212 dígitos :)

Batman e a escalabilidade

Dia desses eu resolvi colocar os meus gibis em ordem (tarefa bem difícil, dada a quantidade deles). Mas logo no começo eu já apanhei pra conseguir ordenar dois deles. Ambos eram Punisher #1, ambos do Garth Ennis, mas qual veio antes?

Normalmente esse tipo de dúvida eu resolvo consultando o Grand Comics Database, um repositório imenso de informações sobre gibis. Porém, após usar um pouco a busca do site, eu notei que ela, hm, não funciona tão bem.

Por exemplo, um dos meus gibis é o Sergio Aragonés' Groo: Mightier than the Sword. Se você buscar no site por "mightier than sword", ele não acha, porque só faz match exato. Se você buscar por "Groo", também não acha. O nome foi indexado como "Groo:", e sem os dois pontos no final ele não acha. Por fim, se você buscar por "Aragones" ele também não acha, porque você não colocou o acento.

A primeira coisa em que pensei foi "poxa, é simples fazer melhor que isso". Mas falar é fácil né? O jeito é arregaçar as mangas e fazer um search que funcione!

Design

Vamos então aos requisitos: eu quero fazer um buscador que me retorne todas as séries que constem no database do GCD. O buscador precisa ser web-based, escalável e de latência baixa.

Para conseguir latência baixa, o jeito mais fácil é fazer alguma coisa que seja AJAX-like. E pra isso, minha ferramenta predileta é o Google Web Toolkit. Eu até gosto de programar em javascript, acho a linguagem legal, mas o problema é que cada browser implementa javascript de um jeito, e na prática você precisa perder um tempão com as peculiaridades de cada browser.

Com o GWT você escreve a aplicação uma vez só, em Java, e ele então compila seu código, gerando um javascript customizado para cada browser. Quando você faz o request, ele faz o download só do javascript correspondente, assim você não precisa gastar banda baixando código do Firefox se está rodando o Chrome. Daí em em diante é tudo AJAX, ou melhor, AJAJ (por default ele usa JSON ao invés de XML).

Para que a aplicação fosse escalável, eu usei o Google App Engine. Com ele, você pode rodar seu servidor direto na infraestrutura do Google, com todas as vantagens que isso implica. Por exemplo, você não precisa se preocupar em ficar dimensionando quantas máquinas são necessárias, ele faz isso sozinho pra você, adicionando máquinas conforme a carga no servidor começa a subir.

O App Engine ainda tem outra vantagem: ele é gratuito! Quer dizer, é gratuito até um certo limite, mas um buscador de gibis é uma aplicação de nicho, e na prática eu não vou ter tráfego suficiente para precisar passar desses limites. Por exemplo, o limite do database gratuito é 1GB, mais que suficiente pra guardar todos os dados que vão ser buscados.

Implementação

Definida a plataforma, vamos à implementação. O jeito mais fácil de fazer um buscador é usando uma lista invertida. Por exemplo, suponha que o database contém os gibis abaixo:

1 - Green Arrow
2 - Green Lantern
3 - Iron Lantern
4 - Green Lantern: Rebirth

A lista invertida correspondente fica assim:

arrow: 1

green: 1,2,4

iron: 3

lantern: 2,3,4
rebirth: 4

Tendo a lista invertida é fácil fazer a busca: basta procurar na lista invertida cada palavra da minha query, e fazer a intersecção dos conjuntos resultantes. Por exemplo, "Green Lantern" retorna {1,2,4} e {2,3,4}, cuja intersecção é {2,4}, que são os dois gibis que vamos retornar. Se as listas estiverem ordenadas, você pode calcular essa intersecção em O(n), e se não estiverem ordenadas, você ainda pode fazer em O(n) usando um pouco de espaço extra (como?)

Dentro do App Engine, as listas invertidas vão ficar guardadas no Datastore, que é um database escalável construído em cima da BigTable. Vale notar que o Datastore não é um database convencional: ao invés da analogia comum de uma tabela com linhas e colunas, o mais apropriado é pensar nele como se fosse um dicionário de dicionários: cada "linha" agora é a chave do dicionário mais externo, e as "colunas" são as chaves dos dicionários internos. Dessa maneira, não é obrigatório que cada "linha" tenha sempre as mesmas "colunas".

Adicionalmente à lista invertida, nós vamos guardar também uma lista direta. Assim, quando eu descobrir que os gibis a mostrar são o 2 e o 4, eu posso ir na lista direta e pegar mais dados sobre eles, como nome, editora, ano de lançamento, etc.

Depois disso, nós ainda precisamos fazer alguma espécie de ranking pra descobrir em qual ordem vamos mostrar os resultados. Em geral essa é uma das partes mais difíceis de uma search engine, mas eu resolvi simplificar e ordenar por número de edições (então Action Comics, que tem quase 900 edições, vai aparecer antes de Slave Girl Comics, que tem só duas).

Por fim, todo esse trabalho de lista invertida, intersecção, lista direta e ranking pode eventualmente ficar meio pesado, então nós vamos colocar um cache no caminho para guardar todas as queries já feitas, sem precisar repetir trabalho. O App Engine fornece um módulo de Memcache que é ideal para isso.

Preparando os dados

Antes de começar a usar o buscador, precisamos gerar as listas e colocá-las dentro do Datastore. O GCD fornece um arquivão com o dump do MySQL deles (que na verdade é um arquivo texto com 330MB de INSERTs). Com os dados em mãos, eu escrevi um script python que lê esses dados, cria a lista invertida e a lista direta, e faz o bulk upload desses dados para o Datastore.

O único cuidado extra nessa etapa é tomar cuidado com o unicode, já que eu quero que buscas por "mônica" e "monica" retornem o mesmo resultado. Uma maneira seria tirando os acentos no cliente, mas pra que fazer isso online se você pode fazer offline né? Mais fácil, na hora de criar a lista invertida, gerar duas entradas iguais, uma com "mônica" e outra com "monica".

Outra idéia que tive foi retornar as capas dos gibis, e não apenas os nomes deles. Nesse caso o processo foi um pouco mais complicado, porque o GCD não fornece um arquivão com todas as capas, assim como eles fazem com o dump do MySQL. Eles acham que fornecer esse arquivão para download poderia ser uma violação do fair use, então a recomendação é "se você quer as capas, faça um crawling você mesmo".

Bora fazer um crawler então! Usando várias threads em python foi tranquilo, eu baixei os 500MB de capas em poucas horas. Como o Datastore não foi feito para guardar imagens, eu optei por hospedar as capas no meu próprio servidor (www.ricbit.com).

O cliente

Escrever o cliente usando o GWT foi a parte mais difícil. Não que programar em GWT seja difícil, pelo contrário, usando UiBinder é bastante fácil. O problema mesmo é que eu manjo de CSS menos do que gostaria, então demorei um tempão até conseguir o layout desejado. Mas, como esperado, o resultado funciona em tudo quanto é browser, até mesmo no IE!

E para dar um toque especial no layout, eu pedi para a minha esposa fazer umas ilustrações para o logo e para a tela de loading do site, que ficaram bem bacanas (se você precisa de ilustrações para o seu site, fale com a Ila Fox!)

Por fim, eu ainda precisava de um jeito de fazer profiling do site, então fiz com que o servidor me retornasse o tempo gasto em cada etapa do processamento, e usei o Google Chart Tools para apresentar esses dados de uma maneira fácil de interpretar.

Segurança

Você não pode fazer uma aplicação web-based sem pensar na segurança, certo? Mas, felizmente, esse buscador não tem muitos pontos vulneráveis. O único ponto onde o usuário malicioso pode tentar fazer alguma coisa é no campo da query, mas eu não preciso nem sequer sanitizar, porque tudo que a aplicação retorna são os nomes dos gibis que já estão no datastore.

Na verdade, só de rodar a aplicação no App Engine você ainda ganha na faixa um monte de vantagens, entre elas detecção de ataques de DoS. Ele sabe sozinho filtrar ips que estão te bombardeando sem nem precisar configurar nada.

Juntando as peças

No final, o sistema completo ficou assim:

Se você quiser testar a aplicação, basta clicar no screenshot abaixo:

Agora finalmente dá para achar aquele gibi do Groo :)

Profiling

Com o sistema completo e funcionando já dá pra fazer o profiling. Para isso, basta digitar "debug:" antes de uma query qualquer. Por exemplo, se você fizer "debug: Groo" logo após o deploy do servidor, o resultado vai ser esse:

O App Engine faz resource management de maneira bem agressiva, se você não tiver nenhuma query em 10min, ele derruba sua instância. Daí, a primeira query que você fizer tem o overhead de levantar o serviço novamente, que nesse caso foi de quase 5 segundos. Mas se você fizer em seguida a mesma query, ela vai ser bem mais rápida:

Bem melhor né? Alguns poucos milisegundos no memcache, o resto é só o tempo de trafegar o resultado do servidor pro seu cliente. No total, pouco menos de 350ms de roundtrip, o que é bastante rápido.

Mas essa é uma query fácil. Uma das queries mais difíceis é "debug: Batman", porque existem muitos, muitos gibis do Batman. Vejamos:

O tempo é praticamente todo gasto baixando os dados dos gibis da lista direta, quase 3 segundos. Mesmo considerando que é um batch request de 824 gibis, três segundos é meio lento. Mas ao usar o Datastore, o importante é notar que esse tempo é independente da carga. Você até consegue uma solução customizada mais rápida que isso, mas sem tomar muito cuidado, o tempo de cada request vai aumentar se o número de usuários simultâneos crescer. Com o Datastore, o tempo é sempre o mesmo, e você está imune a baleiadas.

Naturalmente, isso só vale para a primeira vez. Se outro usuário buscar por Batman logo em seguida, ele vai pegar o resultado do memcache:

Agora deu só 500ms, e o tempo todo é praticamente só o overhead do download dos 824 gibis.

Outra query curiosa é "debug: Captain America Comics". Existem inúmeros gibis com captain (Captain America, Captain Marvel, etc), inúmeros com america (Captain America, Justice League of America, etc), e um outro sem número de Comics, então o tempo agora não é dominado pela lista direta:

Dessa vez dá pra ver o tempo que ele gasta fazendo a intersecção das listas invertidas (e note que até agora nós nem vimos tempo algum com ranking, é desprezível). O tempo total foi de 800ms, mesmo sem o memcache.

O source

É claro que isso tudo foi uma visão bem high-level do Gibit, tem um monte de detalhes curiosos que eu omiti para não deixar o post muito grande. Se você quiser ver a implementação, o código é open source e está disponível no Google Code:

Source code do Gibit

Se você tiver alguma dúvida específica, é só deixá-la nos comentários que eu tento responder. Enjoy :)