segunda-feira, 31 de janeiro de 2022

Automatizando código para salvar arquivos

Quando estamos realizando alguma tarefa, normalmente iterativa, às vezes queremos salvar nosso progresso atual antes de toda a tarefa terminar. Isso evita que interrupções no software façam com que nossos resultados sejam perdidos. Por exemplo, você está realizando um Gibbs sampling e não quer perder suas estimativas caso algo aconteça antes do algoritmo terminar. Sabemos que modelos complexos em Inferência Bayesiana necessitam de mais tempo para realizar inferência por amostragem.
Suponha que você colocou o algoritmo para rodar 100 mil iterações. Mas, a cada 10 mil, quer salvar seu progresso, para evitar perda de resultados. Como fazemos isso?

Primeiro, vamos aprender a salvar um arquivo. Aqui, vou ensinar a usar a função saveRDS(). Essa função salva um objeto do R em um arquivo com extensão .rds e nome especificado pelo usuário. Vejamos um exemplo:
dados = c(1,20,12,90,5);
saveRDS(dados,"teste01.rds");

dados recebe um vetor com os números acima. Em seguida, salvamos o que tem na variável dados no arquivo "teste01.rds". Dessa forma, no diretório atual (que no caso, é o diretório do meu projeto) o arquivo é criado. É importante não esquecer de colocar a extensão do arquivo.

Agora, vamos ler o arquivo e ver o que tem nele.

dadosr = readRDS("teste01.rds")

> dadosr
[1]  1 20 12 90  5

Como esperado, dadosr recebe o que tem no arquivo "teste01.rds", que é o vetor de dados que salvamos inicialmente. Mas, apenas vetores podem ser salvos dessa forma? A resposta é não! Qualquer objeto pode ser gravado dessa maneira. Vejamos um exemplo com lista:

dados_aluno = list(nome = "João", disciplina = "Matematica", nota = 8.7)
saveRDS(dados_aluno,"teste02.rds");

Pronto, a lista criada que armazena o nome, disciplina e nota agora está salva no arquivo "teste02.rds". É claro que há maneiras mais eficientes para armazenar esse tipo de dado, mas aqui é apenas uma ilustração. Se lemos esse arquivo, obtemos:

dados_aluno_r = readRDS("teste02.rds")

> dados_aluno_r
$nome
[1] "João"

$disciplina
[1] "Matematica"

$nota
[1] 8.7

Agora que você sabe salvar um arquivo, é trivial automatizar gravação de resultados. Normalmente o fazemos em algoritmos iterativos e adicionamos o comando de gravação dentro da estrutura de repetição. Obviamente, temos que ter um controle de quanto em quanto tempo queremos gravar os resultados, pois, se gravamos a cada iteração, o algoritmo pode ficar extremamente lento. 

No exemplo mencionado, podemos optar por salvar um arquivo dos resultados a cada 10 mil iterações, até que as 100 mil sejam finalizadas. Ou podemos optar por gravar a cada 1 mil iterações. Esse "tempo" depende do problema em mãos. Se é algo que demora bastante (modelos muito complexos), salvar a cada poucas iterações é interessante. Agora, se o algoritmo é relativamente rápido, podemos dar um salto maior para a gravação. É claro que, se o algoritmo é rápido (algo factível sem o risco de perda de resultados), não precisamos salvar resultados preliminares.

Nesse post, faremos um exemplo banal para ilustrar a importância da automatização de código para gravar arquivos. Vamos gerar duas variáveis aleatórias e o resultado será gravado a cada mil iterações de um total de 10 mil. Vejamos o código abaixo:

result = matrix(ncol=2);
result = result[-1,];

for(i in 1:10000){
  x1 = rnorm(1); # gera uma observação da normal padrão
  x2 = runif(1); # gera uma observação da Uniforme(0,1)
  
  result = rbind(result,cbind(x1,x2));
  
  if(i%%1000 == 0) saveRDS(result,paste0("save_files/result_temp",i,".rds"));
}
saveRDS(result,"result.rds");

Primeiro, criamos a variável result, que armazenará uma matriz. Como ela é iniciada com uma linha contendo NA's, excluímos essa linha. Em seguida, iniciamos uma iteração de tamanho 10 mil. Para cada iteração, geramos um valor de uma distribuição normal padrão e um valor de uma uniforme(0,1) e esse resultado é adicionado a uma linha na matriz. Agora a parte interessante:

if(i%%1000 == 0) saveRDS(result,paste0("save_files/result_temp",i,".rds"));

O que essa linha de código está fazendo? Verificamos se passaram-se mil iterações. Caso positivo, salvamos nosso resultado prévio. Note que estou salvando o arquivo na pasta "save_files", com um nome que depende do valor de i. Ou seja, a cada mil iterações, vamos salvar nosso resultado com o nome "result_temp1000.rds", "result_temp2000.rds" e assim por diante. Se executarmos o código acima, veremos no nosso diretório os seguintes arquivos:

Esses arquivos foram salvos para você visualizar que o processo está funcionando. Se modificarmos o código para salvar o arquivo com o mesmo nome, a cada gravação o arquivo será sobrescrito. Isso é útil, já que queremos sempre o arquivo mais atualizado. Então, modificamos o código para:

if(i%%1000 == 0) saveRDS(result,"save_files/result_temp.rds");

Agora, os resultados prévios serão sempre salvos no arquivo "result_temp.rds". Eu particularmente gosto de colocar a palavra "temp" no nome para indicar que é um arquivo temporário ou de resultados prévios. Após a execução de todo o algoritmo, salvamos o resultado no arquivo "result.rds", que contém o resultado completo.

Note que, após finalizar o for, o arquivo "result_temp.rds" também armazenará o resultado completo. Normalmente isso depende se a última iteração do algoritmo vai salvar ou não os resultados prévios e depende da estrutura do nosso algoritmo. Varia de problema para problema. De qualquer modo, analise o seu problema e veja se é necessário salvar um arquivo após a execução da estrutura de repetição. No nosso caso, não é necessário e seria mal uso de memória física e de tempo, principalmente se o arquivo for muito grande. Sempre pense nesses detalhes!

Espero que tenha gostado da aula.

Até a próxima aula!

sexta-feira, 28 de janeiro de 2022

Automatizando código para criação de gráficos

Suponha que você tenha um data frame (ou uma matriz) com notas de alunos, em que cada coluna representa uma matéria. Suponha também que você queira fazer um gráfico (dispersão, histograma,...) para as notas dos alunos de cada matéria. Se há dez disciplinas, você precisa de dez blocos de código para construir esses gráficos? A resposta é não!

Sabendo identificar os padrões dos gráficos que você quer criar, automatizar a criação destes poupa (e muito) seu trabalho. Hoje, vamos aprender a construir gráficos que são semelhantes, sem a necessidade de criar um bloco de código para cada gráfico.

Voltando ao exemplo, vamos construir um data frame com notas de 100 alunos em 5 disciplinas.

dados = matrix(data = ceiling(runif(500)*100), nrow=100,ncol=5);
colnames(dados) =  c("Matematica","Fisica","Biologia","Quimica","Portugues");
dados = as.data.frame(dados);

No caso, eu gerei as notas dos alunos através de uma distribuição Uniforme, para termos o exemplo. Podemos visualizar o que temos no data frame:

head(dados)

> head(dados)
  Matematica Fisica Biologia Quimica Portugues
1         26      5      100      35        46
2         99     61       56      58        14
3         20     93       50      51        56
4         50     57       73      38        27
5         67     77       21      89        41
6         56     30       41       1        38

Como não fixei uma semente, os valores para você estarão diferentes. De fato, são irrelevantes para o propósito desse post, então não se preocupe com isso.

Podemos fazer um gráfico com todas as variáveis (que estão nas colunas), digitando:

plot(dados)

O que nos dá o seguinte gráfico:

Esse gráfico é excelente para uma análise exploratória se quisermos visualizar como as variáveis se relacionam duas a duas. Mas não é nosso caso. Vamos construir um histograma para cada disciplina. Por exemplo, o código abaixo nos dará o seguinte gráfico:

hist(dados$Matematica,xlab='Nota',ylab='Frequência',main='Matematica')
Até aqui nada de novo. O que queremos é construir o histograma acima para cada disciplina sem a necessidade de escrever a linha de código para cada uma. Isso é possível, como eu disse, identificando os padrões dos gráficos que você quer criar.

Por exemplo, para construir o histograma da disciplina de "Fisica", basta alterar dados$Matematica para dados$Fisica (e os respectivos labels no gráfico). Porém, dessa forma, fica complicado automatizar nosso código. Devemos escrever "Fisica" ou "Biologia" para acessar tais colunas? Não! Podemos simplesmente acessar as colunas, isto é, para acessar as notas de matemática, digitamos dados[,1], para as notas de física, dados[,2]. Agora sim podemos acessar as disciplinas presente nos dados percorrendo as colunas do nosso data frame. O código acima será escrito por


hist(dados[,1],xlab='Nota',ylab='Frequência',main='Matematica')

e dará o mesmo resultado. Você deve ter percebido o padrão aqui. Mudamos a disciplina trocando o número da coluna que queremos acessar em dados. Note que também temos que alterar o título do gráfico. Para isso, basta digitarmos colnames(dados)[1], por exemplo, para acessar o nome da coluna "Matematica".

Vamos fazer o código para gerar os histogramas para as cinco disciplinas. Basta colocarmos a linha de código acima dentro de uma estrutura de repetição, alterando a coluna que será acessada do nosso data frame. Vejamos:

par(mfrow=c(2,3))

for(i in 1:ncol(dados)) hist(dados[,i],xlab="Nota",ylab='Frequência',main=colnames(dados)[i])

A primeira linha de código acima apenas diz para o R fazer uma grade de plots com espaço para 6 gráficos (2 linhas e 3 colunas).

A segunda linha de código é nossa automatização da construção dos histogramas. O que fizemos ali? Para cada coluna do data frame, construímos o histograma referente a coluna (disciplina). O gráfico abaixo é o resultado da execução do código acima.


Interessante, não acha? Apenas com uma linha de código, construímos os histogramas para cada uma das disciplinas dos nossos dados. Se tivéssemos 14 disciplinas, apenas uma linha de código seria suficiente para construir um gráfico para cada uma. Esse é o potencial da automatização de código.

Agora, nem sempre queremos apenas plotar e visualizar o gráfico, mas salvá-los! Vamos fazer isso: salvar nossos histogramas de modo automático. Para isso, utilizaremos a função png (que salva os gráficos em formato .png). Vamos ao código:

diretorio = "imagens/";

for(i in 1:ncol(dados)){
  nome = paste0(diretorio,"grafico_notas_",colnames(dados)[i],".png")
  png(nome,width=480,height=480,pointsize=18)
  
  plot(dados[,i],xlab="Aluno",ylab=colnames(dados)[i],pch=19)
  
  dev.off()
}

Agora a explicação:

  1. Para salvar um arquivo, precisamos saber onde salvar. No caso, o R já está no diretório do meu projeto e quero gravar o arquivo na pasta "imagens". Dessa forma, guardo essa informação na variável diretorio.
  2. Como o nome do gráfico muda para cada disciplina, a variável contendo o nome deve estar dentro da estrutura de repetição. No caso, concatenei o diretorio + "grafico_notas" + disciplina + ".png". Nunca se esqueça de colocar a extensão correta do arquivo. Um exemplo de nome é "imagens/grafico_notas_Matematica.png", ou seja, o arquivo "grafico_notas_Matematica.png" será gravado na pasta imagens.
  3. Para começar a gravação, utilizamos a função png, passando como atributo o nome para gravação. Width informa a largura do gráfico, height informa a altura e pointsize o tamanho dos pontos no gráfico.
  4. Em seguida, plotamos o que queremos salvar. Note que, quando iniciamos a função png, o plot será gravado no arquivo e nada aparecerá no ambiente do R.
  5. Para finalizar a gravação, digitamos dev.off().

Se você executar o código acima, verá que, na pasta "imagens", haverá cinco arquivos, que são os histogramas das notas dos alunos em cada disciplina:


Dessa forma, fica bem mais fácil de construir os gráficos para análise. E vimos apenas um exemplo, mas essa ideia serve para qualquer situação em que você tenha que fazer vários gráficos semelhantes. Um exemplo bem comum é a construção da trajetória das cadeias dos parâmetros em inferência Bayesiana. O gráfico é o mesmo para parâmetros diferentes. Automatizar a construção desses gráficos poupa bastante tempo e trabalho.

Espero que tenha gostado da aula.

Até a próxima aula!

quarta-feira, 26 de janeiro de 2022

Conway's Game of life - Parte III

Nos posts anteriores, implementamos o jogo da vida (Conway's Game of life). Você pode ler aqui a Parte I e a Parte II. Um exemplo de reprodução do jogo é mostrado abaixo.


O problema que tínhamos na nossa implementação é ter que rodar a função para cada passo do algoritmo. Como discutido anteriormente, se queremos avançar n passos, poderíamos usar uma estrutura de repetição que executasse nosso algoritmo n vezes. Qual o problema dessa ideia? Devemos saber exatamente quantos passos queremos dar no jogo. Além disso, uma das maravilhas da simulação do jogo da vida é exatamente visualizar as coisas acontecendo! Para não perdermos esse brilho, vamos implementar uma estrutura de repetição em que apertamos Enter cada vez que queremos avançar no algoritmo e digitamos 'p' para parar.

Vejamos como fica nosso código:

para_alg = F;
while(para_alg==F){
  estado_atual = novo_estado;
  vizinhos = calcula_vizinhos(estado_atual);
  novo_estado = proximo_estado(estado_atual,vizinhos);
  mostra_pop(novo_estado);
  
  varc = readline(prompt = "\n Digite 'p' para parar ou Enter para continuar: ");
  if(varc=='p') para_alg=T;
}

Vou explicar o que está acontecendo. Declaramos uma variável flag, que irá indicar quando o comando while deve parar de executar. Enquanto essa variável tiver valor False, nosso bloco de código será executado.

As primeiras três linhas dentro do comando while nós já vimos. De fato, foram implementadas no post anterior. A função mostra_pop é uma função para criar os gráficos do jogo, como os três primeiros gráficos desse post. Deixarei essa função no final desse post, mas desafio você a tentar escrevê-la sem minha ajuda. Será um bom treino.

As próximas duas linhas de código são a parte que queríamos: pergunta ao usuário se quer continuar ou parar. Caso o usuário queira parar, digita 'p' e o comando while para de ser executado. Dessa forma, podemos simplesmente pressionar Enter para continuar a execução do nosso jogo até ficarmos satisfeitos e querermos pará-lo. Note que, na verdade, qualquer texto que digitarmos que não seja 'p' fará o algoritmo dar mais um passo (mas você pode tratar isso, deixo como um treino).

Podemos melhorar isso? Sim. Eu disse no começo do post que é fácil implementar essa dinâmica para n passos, porém perderíamos o brilho de ver o algoritmo executando, já que os comandos são executados muito rapidamente. Porém, podemos utilizar uma função que faz com que o R espere um tempo antes de executar o próximo bloco de código (que, no nosso caso, seria fazer outra iteração). Utilizaremos a função Sys.sleep(), em que passamos como argumento o tempo (em segundos) em que o R ficará em suspensão. Vejamos como fica nosso código:


n=10;
for(i in 1:n){
  estado_atual = novo_estado;
  vizinhos = calcula_vizinhos(estado_atual)
  novo_estado = proximo_estado(estado_atual,vizinhos);
  mostra_pop(novo_estado)
  
  Sys.sleep(1);
}

Nesse caso, veremos visualmente o algoritmo executar os 10 passos, tomando 1 segundo a cada passo. É o suficiente para vermos de forma clara a simulação do jogo acontecer. O problema dessa implementação é que, se quisermos parar em um determinado passo, não conseguimos. Paramos apenas após os n passos que definimos inicialmente.

Agora sim, temos maneiras fáceis de visualizar a simulação do jogo da vida.

Falta mostrar o código da função mostra_pop. E aí, você conseguiu construir o gráfico do jogo? Vai aí a minha resolução:

mostra_pop = function(estado,titulo=NULL){
  nrow_ = nrow(estado);
  ncol_ = ncol(estado);
  
  plot(0,0,type='n',ylim=c(nrow_+0.5,0.5),xlim=c(0.5,ncol_+0.5),
       xlab='',ylab='',main=titulo,xaxs='i',yaxs='i');
  abline(v=seq(-0.5,(ncol_+0.5),by=1),
         h=seq(-0.5,(nrow_+0.5),by=1));
  
  id_1 = which(estado==1);
  linha = id_1%%nrow_;
  coluna = ceiling(id_1/nrow_);
  
  rect(coluna-0.5,linha-0.5,coluna+0.5,linha+0.5,col=1);
}

Basicamente o código se baseia nas funções plot e rect. Como já falamos delas por aqui, não há necessidade de explicá-las. Vou focar nas seguintes três linhas de código:


  id_1 = which(estado==1);
  linha = id_1%%nrow_;
  coluna = ceiling(id_1/nrow_);

A primeira linha vai identificar as células que estão vivas na gride. O cuidado aqui é que estado se refere a uma matriz e a função which irá retornar os índices como se a matriz fosse um vetor. Por exemplo, em uma matriz 2x2, temos quatro células, e a função which irá retornar os índices (se existirem na condição) de 1 a 4. Dessa forma, precisamos identificar na matriz quais são esses índices.

Para isso, temos que saber como a função which identifica tais índices. Ela faz isso 'transformando' a matriz em um vetor, em que as colunas são concatenadas. No exemplo da matriz 2x2, o vetor ficaria da seguinte forma: ([1,1],[2,1],[1,2],[2,2]) que corresponde aos índices (1,2,3,4). 

Sabendo disso, para identificar a que linha determinado índice se refere, basta tomarmos o resto da divisão desse índice pelo número de linhas da matriz. Por exemplo, o índice 3 (que é a célula [1,2]) se refere a linha 1, que é o resto da divisão de 3 por 2.

Para identificar a coluna, basta arredondar para cima a divisão do índice pelo número de linhas. No caso, 3 dividido por 2 é 1,5 e, arredondando para cima, identificamos a segunda coluna. A ideia é como se tivéssemos percorrido uma coluna e meia, o que indica que estamos parados na segunda coluna. Dessa forma, identificamos, na matriz, quais células estão vivas e desenhamos o retângulo preenchido no gráfico.

Espero que tenham gostado.
Até a próxima aula!

segunda-feira, 24 de janeiro de 2022

Conway's Game of life - Parte II

No post anterior, começamos a programar o jogo da vida (Conway's Game of life). Você pode ler aqui a Parte I. Vamos dar continuidade à nossa implementação.

Já construímos a função que calcula o número de vizinhos de cada célula no estado atual. Dados o estado atual e o número de vizinhos, precisamos aplicar as regras do jogo. Vamos relembrá-las:

  1. Qualquer célula viva com menos de dois vizinhos vivos morre de solidão.
  2. Qualquer célula viva com mais de três vizinhos vivos morre de superpopulação.
  3. Qualquer célula morta com exatamente três vizinhos vivos se torna uma célula viva.
  4. Qualquer célula viva com dois ou três vizinhos vivos continua no mesmo estado para a próxima geração.

Vamos criar uma nova matriz, chamada de novo_estado, em que representa a configuração do jogo no próximo passo. Inicialmente, essa matriz será igual a matriz do estado atual:

  novo_estado = estado_atual;

Agora, para cada célula, as regras de (1) a (4) acima citadas devem ser aplicadas. Primeiro, verificamos se a célula está viva ou morta:

  if(estado_atual[i,j]==0){ #celula morta
    ...
  }else{ #celula viva
    ...
  }

Se a célula está morta, aplicamos a regra (3), caso contrário, aplicamos as restantes. Dessa forma, temos o seguinte código:

  if(estado_atual[i,j]==0){ #celula morta
    if(vizinhos[i,j]==3) novo_estado[i,j] = 1;
  }else{ #celula viva
    if(vizinhos[i,j] < 2 || vizinhos[i,j] > 3) novo_estado[i,j] = 0;
  }

Note que não precisamos implementar os casos em que a célula está morta e continua morta ou está viva e continua viva. Isso porque a matriz de novo estado é inicializada com os valores da matriz de estado atual. Se iniciássemos a matriz de novo estado completa de zeros, por exemplo, precisaríamos fazer essas trocas e, consequentemente, nosso algoritmo perderia eficiência.

Podemos colocar nosso código dentro de uma função (o que é recomendável):

  proximo_estado = function(estado_atual,vizinhos){
    ncol_ = ncol(estado_atual);
    nrow_ = nrow(estado_atual);
    novo_estado = estado_atual;

    for(i in 1:nrow_){
      for(j in 1:ncol_){
        if(estado_atual[i,j]==0){ #celula morta
          if(vizinhos[i,j]==3) novo_estado[i,j] = 1;
        }else{
          if(vizinhos[i,j]<2 || vizinhos[i,j]>3) novo_estado[i,j] = 0;
        }
      }
    }
    return(novo_estado);
  }

Nosso jogo, tecnicamente, está pronto. Demos um passo inicial para o algoritmo, calculamos os vizinhos e implementamos o novo estado. Vamos verificar se funciona? Vamos utilizar uma gride de 6 por 6.

ncol_ = 6;
nrow_ = 6;
estado_atual = matrix(0,ncol=ncol_,nrow=nrow_)

pop2 = expand.grid(1:3,2:5)
estado_atual[pop2[,1],pop2[,2]] = 1
vizinhos = calcula_vizinhos(estado_atual)
novo_estado = proximo_estado(estado_atual,vizinhos);

> estado_atual
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    1    1    1    1    0
[2,]    0    1    1    1    1    0
[3,]    0    1    1    1    1    0
[4,]    0    0    0    0    0    0
[5,]    0    0    0    0    0    0
[6,]    0    0    0    0    0    0
> vizinhos
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    2    3    5    5    3    2
[2,]    3    5    8    8    5    3
[3,]    2    3    5    5    3    2
[4,]    1    2    3    3    2    1
[5,]    0    0    0    0    0    0
[6,]    0    0    0    0    0    0
> novo_estado
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    1    0    0    1    0
[2,]    1    0    0    0    0    1
[3,]    0    1    0    0    1    0
[4,]    0    0    1    1    0    0
[5,]    0    0    0    0    0    0
[6,]    0    0    0    0    0    0

Na figura abaixo, vemos o estado atual (esquerda) e o novo estado (direita).

Para executar o próximo passo do jogo da vida, o novo estado vira o estado atual e procedemos da mesma forma, isto é:
estado_atual = novo_estado;
vizinhos = calcula_vizinhos(estado_atual);
novo_estado = proximo_estado(estado_atual,vizinhos);

A figura abaixo representa os próximos dois passos do algoritmo:


Nosso código está funcionando! O próximo passo é implementar uma forma de o algoritmo fazer alguns passos sem precisarmos rodar o bloco de código manualmente a cada passo. Para n passos, é fácil: basta colocar o código em uma estrutura de repetição que irá executar n vezes. Mas e se quisermos ter controle da execução, por exemplo, apertar Enter para o algoritmo executar até ficarmos satisfeitos e depois parar a execução? Esse problema é parecido com o do post de preencher retângulos no plot (que você pode ver aqui), mas vamos abordá-lo no próximo post.

Até a próxima aula!

sexta-feira, 21 de janeiro de 2022

Conway's Game of life - Parte I

O jogo da vida é um autômato celular criado em 1970 pelo matemático John Horton Conway. O jogo simula, através de regras simples, comportamentos em grupos de seres vivos. A cada passo do algoritmo, uma determinada célula pode viver, morrer ou continuar no seu estado atual. Um exemplo do jogo da vida é mostrado abaixo:
Imagem retirada da página do jogo da vida na wikipédia.

As regras do jogo da vida são as seguintes:
  1. Qualquer célula viva com menos de dois vizinhos vivos morre de solidão.
  2. Qualquer célula viva com mais de três vizinhos vivos morre de superpopulação.
  3. Qualquer célula morta com exatamente três vizinhos vivos se torna uma célula viva.
  4. Qualquer célula viva com dois ou três vizinhos vivos continua no mesmo estado para a próxima geração.

No post de hoje, vamos reproduzir o jogo da vida no R. Verificando as regras acima, percebemos que precisamos de três coisas: o estado atual do algoritmo, o número de vizinhos para cada célula do estado atual, e o estado no próximo passo do algoritmo (novo estado). Assim que o novo estado é computado, ele passa a ser o estado atual e o algoritmo procede da mesma forma.

Para implementar o jogo da vida, vamos definir nossa gride de células. 
ncol_ = 11;
nrow_ = 11;

estado_atual = matrix(0,ncol=ncol_,nrow=nrow_)

A matriz estado_atual define qual o estado atual da simulação. No momento, não temos nenhuma célula viva. Temos que definir um valor inicial para o algoritmo.

pop1 = expand.grid(2:6,5:10)
estado_atual[pop1[,1],pop1[,2]] = 1

Dessa forma, informamos que as células da gride definidas em pop1 estão vivas. A Figura abaixo apresenta o estado inicial do nosso algoritmo, em que as células preenchidas representam as células vivas.

Não se assuste ao perceber que a visualização gráfica está diferente da matriz:
> estado_atual
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]    0    0    0    0    0    0    0    0    0     0     0
 [2,]    0    0    0    0    1    1    1    1    1     1     0
 [3,]    0    0    0    0    1    1    1    1    1     1     0
 [4,]    0    0    0    0    1    1    1    1    1     1     0
 [5,]    0    0    0    0    1    1    1    1    1     1     0
 [6,]    0    0    0    0    1    1    1    1    1     1     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     0

Isso ocorre porque eu defini que as linhas estariam representadas no eixo x e as colunas, no eixo y. É claro que você pode mudar isso, dependendo da visualização gráfica que quer. Para o algoritmo em si, não faz diferença, já que trabalhamos com as matrizes.

Dada a matriz com o estado inicial, precisamos calcular os vizinhos (vivos) de cada célula. Para isso, vamos contar os vizinhos de cada célula e armazenar esse informação em uma matriz de vizinhos, em que vizinhos[1,1] terá como valor o número de vizinhos da célula [1,1] da nossa gride. Vamos criar nossa matriz de vizinhos:

vizinhos = matrix(0,ncol=ncol_,nrow=nrow_)

Note que as algumas células possuem número de células adjacentes diferentes de outras. Por exemplo, a primeira célula da matriz [1,1] possui como células adjacentes as células [1,2], [2,1] e [2,2]. Já a célula [2,2] possui oito células adjacentes. Logo, a implementação da contagem de vizinhos vivos deve levar isso em conta. Então, temos que fazer essa computação com alguns cuidados:

  1. as células dos extremos da matriz possuem apenas 3 células adjacentes,
  2. as células das bordas da matriz (exceto os extremos) possuem 5 células adjacentes,
  3. as células do interior da matriz possuem 8 células adjacentes.

Para contar quantos vizinhos vivos cada célula possui, vamos utilizar uma estrutura de repetição. Utilizaremos um for duplo, ou seja, um dentro do outro. O primeiro percorrerá as linhas da matriz e o segundo, as colunas:

for(i in 1:nrow_){
  for(j in 1:ncol_){
    ...
  }
}

Logo, i indica a linha atual e j indica a coluna atual. Usaremos essas variáveis para acessar a matriz de estado atual e computar os vizinhos vivos da célula [i,j]

Para as células do interior da matriz, as células adjacentes seguem o mesmo padrão: linha i-1, linha i+1, coluna j-1, coluna j+1. Vamos implementar:

for(i in 1:nrow_){
  for(j in 1:ncol_){
    aux = as.matrix(expand.grid((i-1):(i+1),(j-1):(j+1)));
    vizinhos[i,j] = sum(estado_atual[aux])-estado_atual[i,j];
  }
}

Perceba que pegamos um bloco de nove células: as oito células adjacentes à célula [i,j] e a própria célula [i,j]. Dessa forma, quando somamos os valores do estado atual desse bloco, precisamos subtrair o valor da célula [i,j] para obtermos apenas os vizinhos vivos. Parece fácil, correto? Para as células do interior da matriz sim, mas para as células na borda da matriz devemos ter cuidado extra.

Por exemplo, as quatro células dos extremos da matriz precisam de um código diferente para o cálculo de vizinhos, apesar de possuírem a mesma quantidade de células adjacentes. Se pegamos a célula [1,1], as células adjacentes estão na próxima coluna e na próxima linha. Agora, se pegamos a célula [nrow_,1], suas células adjacentes estão na próxima coluna e na linha anterior. Percebem o cuidado que devemos ter aqui? E isso se aplica à todas as células da borda da matriz!

Vamos passo a passo. Como o for da linha é o externo, vamos começar identificando se a linha atual está no interior da célula ou pertence à uma das bordas:

if(i>1 && i<nrow_){
  ...
}else if(i==1){ #primeira linha
  ...
}else if(i==nrow_){ #ultima linha
  ...
}

Agora, para cada situação da linha (ou seja, dentro de cada if do código acima), devemos verificar se a coluna atual está nas bordas ou não, ou seja, se j = 1 ou ncol_ ou está entre esses valores.

if(j>1 && j<ncol_){
   ...
}else if(j==1){ #primeira coluna
   ...
}else{#ultima coluna
   ...
}

Agora, o que colocar dentro de cada condição? Se a célula está no interior da matriz, nós já fizemos o código anteriormente. Mas, se está nas bordas, devemos computar os vizinhos de maneiras particulares. Nesse momento, você deve ter entendido como calcular esses vizinhos para essas células. Dessa forma, abaixo coloco a função completa, de fácil entendimento.

calcula_vizinhos = function(estado_atual){
  ncol_ = ncol(estado_atual);
  nrow_ = nrow(estado_atual);
  vizinhos = matrix(0,ncol=ncol_,nrow=nrow_)
  
  for(i in 1:nrow_){
    for(j in 1:ncol_){
      if(i>1 && i<nrow_){
          if(j>1 && j<ncol_){
              aux = as.matrix(expand.grid((i-1):(i+1),(j-1):(j+1)))
          }else if(j==1){ #primeira coluna
              aux = as.matrix(expand.grid((i-1):(i+1),(j):(j+1)))
          }else{#ultima coluna
              aux = as.matrix(expand.grid((i-1):(i+1),(j-1):(j)))
          }
      }else if(i==1){ #primeira linha
          if(j>1 && j<ncol_){
              aux = as.matrix(expand.grid((i):(i+1),(j-1):(j+1)))
          }else if(j==1){#primeira coluna
              aux = as.matrix(expand.grid((i):(i+1),(j):(j+1)))
          }else{#ultima coluna
              aux = as.matrix(expand.grid((i):(i+1),(j-1):(j)))
          }
      }else if(i==nrow_){ #ultima linha
          if(j>1 && j<ncol_){
              aux = as.matrix(expand.grid((i-1):(i),(j-1):(j+1)))
          }else if(j==1){#primeira coluna
              aux = as.matrix(expand.grid((i-1):(i),(j):(j+1)))
          }else{#ultima coluna
              aux = as.matrix(expand.grid((i-1):(i),(j-1):(j)))
          }
      }
      vizinhos[i,j] = sum(estado_atual[aux])-estado_atual[i,j]
    }
  }
  return(vizinhos);
}

Resumindo a função acima: dada uma determinada célula, pegamos um bloco de células, que depende da posição dessa célula na matriz, e calculamos seus vizinhos (somando os valores do bloco e retirando o valor da própria célula). Nossa função, então, retorna a matriz de números de vizinhos do estado atual do algoritmo do jogo da vida.

Vejamos como funciona para nosso exemplo:

> vizinhos = calcula_vizinhos(estado_atual)
> vizinhos
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]    0    0    0    1    2    3    3    3    3     2     1
 [2,]    0    0    0    2    3    5    5    5    5     3     2
 [3,]    0    0    0    3    5    8    8    8    8     5     3
 [4,]    0    0    0    3    5    8    8    8    8     5     3
 [5,]    0    0    0    3    5    8    8    8    8     5     3
 [6,]    0    0    0    2    3    5    5    5    5     3     2
 [7,]    0    0    0    1    2    3    3    3    3     2     1
 [8,]    0    0    0    0    0    0    0    0    0     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     0
> estado_atual
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]    0    0    0    0    0    0    0    0    0     0     0
 [2,]    0    0    0    0    1    1    1    1    1     1     0
 [3,]    0    0    0    0    1    1    1    1    1     1     0
 [4,]    0    0    0    0    1    1    1    1    1     1     0
 [5,]    0    0    0    0    1    1    1    1    1     1     0
 [6,]    0    0    0    0    1    1    1    1    1     1     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     0

Já deu um trabalho até aqui, certo? Continuamos a implementação do algoritmo no próximo post.

Até a próxima aula!

quarta-feira, 19 de janeiro de 2022

Preenchendo retângulos em um plot

Imagine que você tenha uma gride de células em um plot, como mostrado na figura abaixo à esquerda. Agora suponha que você queira selecionar alguns desses retângulos, podendo preenchê-los de outra cor. Um maneira de fazer seria escrever na mão pontos que estão dentro dos retângulos de sua escolha, por exemplo, os pontos (1,5; 2,5) e (8,5;8,5) correspondem aos retângulos preenchidos na figura à direita.





Se você quer selecionar vários retângulos, escrever os pontos à mão é ineficiente (e extremamente cansativo). Por isso, hoje vamos escrever uma função para preencher o retângulo que queremos quando clicamos dentro dele.

Faremos isso utilizando a função locator(). Essa função pega as coordenadas de um clique no gráfico, ou seja, ao clicar em uma área do gráfico, a função retornará as coordenadas do ponto em que houve o clique. Por exemplo:
> locator(1)
$x
[1] 0.3500569

$y
[1] 5.517714

Com as coordenadas em mãos, sabemos que o retângulo escolhido foi o que está entre 0 e 1 no eixo x e entre 5 e 6 no eixo y. Desse modo, sabemos qual retângulo preencher. Como faremos isso? Da seguinte maneira: desenhamos um retângulo (com cor de preenchimento) que tem ponto central igual a 0,5 e 5,5. Obviamente, esse ponto central é particular da gride que eu criei, caso sua gride tenha intervalos diferentes (por exemplo, retângulos com comprimento igual a 2), ajuste esse ponto.

Para computarmos a coordenada do ponto central do retângulo escolhido, faremos o seguinte: 

  • arredondamos a coordenada x obtida para cima e subtraímos metade do comprimento;
  • arredondamos a coordenada y obtida para cima e subtraímos metade da altura.

No caso, o retângulo tem comprimento e altura iguais a 1, ou seja, subtraímos 0,5 de ambas coordenadas arredondadas:

xr = ceiling(p$x) - 0.5;
yr = ceiling(p$y) - 0.5;

Dado esse centro e as dimensões dos retângulos, estamos prontos para desenhar o retângulo preenchido:

rect(xr-0.5,yr-0.5,xr+0.5,yr+0.5,col=1)

Como o código está solto, vamos colocá-lo dentro de uma função. Além de boa prática, nos ajuda a verificar se a função está funcionando corretamente.

desenha = function(){
  p = locator(1)
  xr = ceiling(p$x) - 0.5;
  yr = ceiling(p$y) - 0.5;
  rect(xr-0.5,yr-0.5,xr+0.5,yr+0.5,col=1)
}

Se executarmos a função acima e clicarmos dentro de um retângulo, ele será preenchido com a cor preta (col = 1). Note que, como parâmetro da função locator(), passamos o valor 1. Isso significa que queremos apenas um ponto. Poderíamos passar n como argumento para podermos colorir n pontos. Vamos testar? A única alteração é passar um valor n para a função desenha() que será utilizado na função locator().

desenha = function(n = 1){
  p = locator(n)
  ...
}

No caso, colocamos n = 1 no argumento da função desenha() para atribuirmos um valor padrão, quando o usuário não passa esse argumento. Temos dois "problemas" na função acima:

  1. o preenchimento ocorre apenas depois de selecionarmos todos os n retângulos,
  2. devemos saber exatamente quantos retângulos queremos colorir.

Não seria mais fácil colorir os retângulos, um por um, até ficarmos satisfeitos e depois parar a função? Muito mais prático, não? Vamos implementar! Para isso, voltamos a passar apenas um ponto na função locator().

Vamos pensar... queremos executar a função desenha() completa uma vez, outra vez, e outra vez até termos um sinal de parada. O que isso nos lembra? A estrutura de repetição, é claro. Então, vamos executar a função desenha() até que um critério de parada seja satisfeito. Vamos considerar que a parada ocorre quando a função locator() não retorna coordenadas, ou seja, não há clique de pontos. Isso ocorre quando o usuário aperta ESC enquanto a função locator() está executando.

Agora, se a função desenha() está executando infinitamente, como saber se o usuário terminou sua ação? Simples, quando não há coordenadas para retornar, a função locator() retorna NULL. Dessa forma, verificamos se a variável que armazena o ponto obtido está nula ou não. Vejamos:

desenha = function(){
  p = locator(1);
  if(is.null(p)){
    cat("Seleção de pontos encerrada.\n");
    return(T);
  }else{
    xr = ceiling(p$x) - 0.5;
    yr = ceiling(p$y) - 0.5;
    rect(xr-0.5,yr-0.5,xr+0.5,yr+0.5,col=1)
    return(F);
  }
}

Na função acima, verificamos se p é nulo. Caso positivo, printamos a mensagem "Seleção de pontos encerrada." e a função retorna TRUE. Caso contrário, preenchemos o retângulo escolhido e a função retorna FALSE. Feito isso, basta escrever a estrutura de repetição:

parada_flag = F;
while(parada_flag==F){
  parada_flag = desenha();
}

Utilizamos a variável parada_flag para identificar quando parar a repetição da função desenha(). Quando o usuário pressiona ESC, a função desenha() retorna o valor TRUE e, então, o while() é terminado. Se você testar o código acima, verá que, a cada clique, o retângulo escolhido é preenchido, exatamente do jeito que queríamos. Assim que terminamos, apertamos ESC e a execução é finalizada.

A função está completa, se consideramos o problema proposto inicialmente. Mas e se quisermos guardar a informação de quais retângulos foram preenchidos? Precisamos, então, retornar as coordenadas dos pontos clicados. 

Aqui vai uma observação: se você quer exatamente as coordenadas dos pontos clicados, deve retornar o valor que a variável p armazena. No meu caso, o interesse é apenas identificar os retângulos e, dessa forma, posso armazenar apenas o ponto central deles. É isso o que vamos fazer:

desenha = function(){
  p = locator(1);
  if(is.null(p)){
    cat("Seleção de pontos encerrada.\n");
    p = NULL;
    parada = T;
  }else{
    xr = ceiling(p$x) - 0.5;
    yr = ceiling(p$y) - 0.5;
    rect(xr-0.5,yr-0.5,xr+0.5,yr+0.5,col=1)
    p = c(xr,yr);
    parada = F;
  }
  return(list(p=p,parada=parada));
}

Perceba que agora, além de retornar T ou F, precisamos retornar o ponto clicado. Dessa forma, retornamos uma lista, em que p identifica o ponto central do retângulo ou o valor NULL (caso não haja ponto clicado) e parada retorna a identificação da condição de parada. Como estamos retornando uma lista, devemos modificar nossa estrutura de repetição. Além disso, precisamos armazenar os pontos clicado em algum lugar. Fazemos isso da seguinte forma:

pontos = matrix(ncol=2);
pontos = pontos[-1,]

Ou seja, pontos é uma matriz de duas colunas. Sempre é inicializada com uma linha, caso não haja especificação. Então a segunda linha do código acima exclui a linha da matriz, deixando-a vazia. Assim, ela está pronta para receber as coordenadas dos pontos (primeira coluna para eixo x e segunda coluna para eixo y).

A estrutura de repetição fica:

parada_flag = F;
while(parada_flag==F){
  aux = desenha(1);
  parada_flag = aux$parada;
  pontos = rbind(pontos,aux$p);
}

Perceba que usei uma variável auxiliar aux para receber o que a função desenha() retorna e, a partir dela, atualizei a variável parada_flag e incrementei o ponto central na matriz pontos. Dessa forma, armazenamos as coordenadas dos pontos clicados (no meu caso, do ponto central do retângulo) cada vez que escolhemos um retângulo para preencher.

Vamos testar nosso código? Meu exemplo ficou da seguinte maneira:

> pontos = matrix(ncol=2);
> pontos = pontos[-1,]
> parada_flag = F;
> while(parada_flag==F){
+   aux = desenha();
+   parada_flag = aux$parada;
+   pontos = rbind(pontos,aux$p);
+ }
Seleção de pontos encerrada.
> head(pontos)
     [,1] [,2]
[1,]  3.5  8.5
[2,]  4.5  8.5
[3,]  5.5  8.5
[4,]  2.5  7.5
[5,]  6.5  7.5
[6,]  1.5  6.5


Perceba que a lógica utilizada aqui serve para colorir qualquer figura geométrica, como triângulos e pentágonos, apenas ajustando o código. Claro que, dependendo da complexidade da figura, o ajuste será mais ou menos difícil.

Você deve estar se perguntando: para que serve o código que fizemos? Tenho duas respostas diretas para isso:
  1. aprender a programar. Não é tão simples implementar essas funções como parece e ajuda na prática de programação.
  2. selecionar o estado inicial do Jogo da vida (Conway's Game of Life). Veremos esse algoritmo no próximo post, mas você pode dar uma olhada nele aqui.

Esperam que tenham gostado da aula.
Até a próxima aula!

segunda-feira, 17 de janeiro de 2022

Customizando gráficos com a função plot

A função plot do R é a função base para se fazer gráficos. Há outras maneiras, por exemplo utilizando o pacote ggplot2, mas quero focar no mais simples, que é a base do R. Desse modo, vamos aprender a alterar tamanho de fonte (texto) e cor, além de formatar a fonte. Para isso, vamos utilizar um problema exemplo que está presente no livro "Noções de Probabilidade e Estatística", de Marcos Magalhães e Antonio Lima, que reproduzo a seguir:

O departamento de vendas de certa companhia ofereceu um curso de atualização a seus funcionários e, para estudar a eficácia do curso, resolveu compara a nota de teste no curso com o volume de vendas, em milhares de unidades, nos seis meses seguintes ao curso. Os resultados estão na tabela abaixo:

Notas 8 9 7 8 6 8 5 5 6 7 4 7 3 5 3
Vendas 14 13 12 13 10 12 11 11 10 12 10 13 10 12 11

Agora, queremos relacionar as duas variáveis. Para isso, vamos plotar as Vendas em função das Notas.


notas  = c(8,9,7,8,6,8,5,5,6,7,4,7,3,5,3)
vendas = c(14,13,12,13,10,12,11,11,10,12,10,13,10,12,11)

plot(notas,vendas)
Perceba que nosso gráfico está simples demais e meio feio para ser apresentado para alguém, não acha? Vamos melhorar isso.

Primeiro, vamos adicionar um título e um texto para os eixos (chamaremos de "labels"). Para isso, utilizamos os parâmetros da função plot:

  • main: adiciona um título;
  • xlab: texto do eixo x;
  • ylab: texto do eixo y.

plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas")
Melhor, não acha? Mas podemos melhorar ainda mais. Vamos aumentar o tamanho de fonte do texto dos eixos. Dessa forma, melhoramos a visualização de texto sem ter que aumentar o gráfico.
Para isso, utilizamos o parâmetro cex da função plot:

  • cex.lab: tamanho de fonte do texto dos eixos;
  • cex.axis: tamanho de fonte dos eixos;
  • cex.main: tamanho de fonte do título.

Vejamos como cada um desses atributos altera o nosso gráfico:

plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas",
     cex.lab=1.5)



plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas",
     cex.axis=1.5)



plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas",
     cex.main=1.5)

Compare esses gráficos com o primeiro. Bem melhor a visualização, certo? Além do tamanho, também podemos alterar o estilo da fonte, utilizando o parâmetro font:

  • font.lab: altera o estilo da fonte do texto dos eixos;
  • font.axis: altera o estilo da fonte dos eixos;
  • font.main: altera o estilo da fonte do título.

Esse parâmetro recebe um valor inteiro para identificar qual estilo o usuário quer, sendo:

  1. texto normal,
  2. negrito;
  3. itálico;
  4. negrito e itálico;
  5. fonte de símbolo.

Vejamos como ficaria nosso gráfico alterando esses atributos:


plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas",
     font.lab=2, font.axis=3, font.main=4)

Note que o texto dos eixos está em negrito, os números dos eixos está em itálico e o título está em negrito e itálico. Também perceba que apenas colocando negrito no texto dos eixos temos uma melhor visualização, a leitura fica mais agradável e direta, sem necessidade de aumentar o gráfico. Em outras palavras, o texto fica mais legível.

Também podemos alterar a cor dos textos. Eu particularmente não utilizo muito esse recurso, mas cabe você decidir se vai usar ou não. Para isso, utilizamos o parâmetro col:

  • col.lab: altera a cor do texto dos eixos;
  • col.axis: altera a cor dos números dos eixos;
  • col.main: altera a cor do título.

Vejamos como fica nosso gráfico com essas alterações:

plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas",
     col.lab="blue",col.main="red",col.axis="green")
Obviamente, o cuidado na escolha das cores é fundamental. Se escolhemos uma cor muito clara em um fundo claro, podemos tornar o respectivo texto ilegível.

O último parâmetro que vamos discutir nesse post é o pch. Esse parâmetro indica a forma dos pontos no gráfico. Por exemplo, estamos utilizando um círculo sem preenchimento como ponto, mas podemos utilizar triângulo, cruz ou um círculo preenchido, entre outros. A imagem abaixo identifica quais valores correspondem a quais formatos (imagem retirada do help do R - ?pch).


Agora que vimos várias formas de alterar atributos nos gráficos do R, vamos refazer nosso gráfico? Veja como ficou:
plot(notas,vendas,xlab="Notas de teste",ylab="Volume de vendas",
     main = "Relação vendas e notas",
     cex.lab=1.2, cex.axis=1.1,cex.main=1.3,
     font.lab=2,
     pch=16)

Compare esse gráfico com o primeiro que fizemos. Agora temos um gráfico apresentável, não acha? E apenas utilizando a função base do R.

Até a próxima aula!



sexta-feira, 14 de janeiro de 2022

Decompondo números em fatores primos

Já que falamos sobre números primos nas últimas aulas, hoje vamos implementar uma função para decompor um dado número em fatores primos. Vamos relembrar como fazemos isso. 

Suponha que queremos decompor o número 30 em fatores primos.

  • Dividimos 30 pelo menor número primo: 2. A divisão é exata, guardamos o número 2 e seguimos com o número 30/2 = 15.
  • Dividimos 15 por 2. A divisão não é exata, então pegamos o próximo número primo.
  • Dividimos 15 por 3. A divisão é exata, guardamos o número 3 e seguimos com o número 15/3 = 5.
  • Dividimos 5 por 3. A divisão não é exata, então pegamos o próximo número primo.
  • Dividimos 5 por 5. A divisão é exata, guardamos o número 5 e seguimos com o número 5/5 = 1.
  • Como não há mais o que dividir (chegamos ao número 1), então o resultado é 2*3*5.

Desse modo, encontrar os fatores primos de um número significa dividir esse número sucessivamente pelo números primos em ordem crescente, guardando os que resultam em divisão exata. Note que, antes de passar para o próximo primo, testamos a divisão pelo primo atual novamente. Por exemplo, o número 8 é dividido por 2 três vezes, resultando em 2*2*2. 

Para implementarmos essa função, necessitamos de uma lista de números primos. Se você não tem uma lista guardada, pode criar uma e salvá-la em um arquivo para uso posterior. É o que faremos aqui, utilizando a função do crivo de Eratóstenes que construímos na última aula:


res = crivo_erat2(200)
saveRDS(res,"lista_primos.rds")
Desse modo, criamos uma lista de números primos e armazenamos na variável res. A segunda linha de código salva o que está armazenado em res em um arquivo chamado "lista_primos.rds". O arquivo será salvo no diretório atual do projeto. Caso queira especificar outro local, digite o caminho inteiro caso o destino esteja fora da pasta do projeto ou as subpastas caso esteja dentro da pasta do projeto.

Supostamente já temos uma lista de primos disponível e o passo acima apenas se refere a como salvar uma lista desse tipo. 

De fato, vamos começar nosso algoritmo carregando a lista de primos salva:

lista_primos = readRDS("lista_primos.rds")

> lista_primos
 [1]   2   3   5   7  11  13  17  19  23  29  31  37  41  43  47  53  59  61  67  71
[21]  73  79  83  89  97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173
[41] 179 181 191 193 197 199

Agora que temos a lista de números primos disponível, vamos implementar a função. Como dito, a decomposição em números primos se dá dividindo o número desejado sucessivamente pelos números primos da lista. O código é bem simples de entender:


fatora_em_primos = function(n,lista_p){
  resto = n;
  ind_p = 1;
  lista_f = NULL;
  while(resto>1){
    if(resto%%lista_p[ind_p]==0) {
      lista_f = c(lista_f,lista_p[ind_p]);
      resto = resto/lista_p[ind_p];
    }else{
      ind_p = ind_p +1;
    }
  }
  return(lista_f)
}

  • A variável resto armazena o número escolhido e atualiza quando é dividido por um número primo.
  • A variável ind_p armazena o índice atual da lista de primos.
  • A variável lista_f armazena os fatores primos do número.
  • A divisão é realizada até que resto seja igual a 1.

Vamos testar nossa função:


numero = 90;
result = fatora_em_primos(numero,lista_primos)

> result
[1] 2 3 3 5

Perceba que nossa função apresenta como resultado todos os fatores primos do número, podendo ser repetidos. No exemplo acima, 90 pode ser decomposto em 2*(3^2)*9. Como indicar quantas vezes o número 3 aparece e reduzir o vetor do resultado? Podemos utilizar a função table().


> table(result)
result
2 3 5 
1 2 1 

Agora sim obtemos uma lista de fatores primos únicos e quantas vezes cada um aparece. Matematicamente falando, os números da linha debaixo do resultado acima correspondem ao expoente do número primo correspondente da fatoração.

Se você não percebeu, há um problema na nossa função: e se a lista de números primos que é fornecida não é suficiente? Isso não é difícil de acontecer, basta pegarmos um número primo maior que o maior número primo da lista de primos. Por exemplo, o número 211 é primo e se executamos nossa função com a lista previamente carregada, obtemos:


> numero = 211;
> result = fatora_em_primos(numero,lista_primos)
Error in if (resto%%lista_p[ind_p] == 0) { : 
  missing value where TRUE/FALSE needed

O erro ocorre porque tentamos acessar um índice depois do fim do vetor de lista_p. Precisamos tratar esse erro.

Para isso, basta adicionar uma condição dentro do comando de repetição while, indicando que se o vetor lista_p foi totalmente percorrido e não obtemos o resultado, então devemos parar.


  while(resto>1){
    if(ind_p > length(lista_p)){
      cat("Lista de primos insuficiente!\n");
      return();
    }
    ...
  }

No código acima, além de pararmos a função e não retornarmos nada, printamos na tela uma mensagem para o usuário: "Lista de primos insuficiente!". Dessa forma, o usuário deve fornecer uma lista de números primos maior.

Se executamos nossa função com o número 211 novamente, obtemos:


> numero = 211;
> result = fatora_em_primos(numero,lista_primos)
Lista de primos insuficiente!

Pronto, agora temos uma função que decompõe um número escolhido em números primos e, caso não haja o fator na lista fornecida, outra lista é requerida.

Até a próxima aula!

quarta-feira, 12 de janeiro de 2022

O crivo de Eratóstenes

Segundo  Coutinho (2005): "o crivo de Eratóstenes é o mais antigo dos métodos para achar primos, e não envolve nenhuma fórmula explícita." Crivo significa peneira e o método consiste em "passar" os números por uma peneira, sobrando apenas os números primos. Dado um número inteiro n positivo, o crivo de Eratóstenes encontra todos os números primos menores que n. Vamos entender como esse método funciona. (lembrando que qualquer número par maior que 2 é composto).

Listamos os números ímpares de 3 a n. Começamos então a riscar os números múltiplos de 3, ou seja, a partir do número 3 riscamos números de 3 em 3. Fazemos o mesmo passo para o número 5, ou seja, riscamos números de 5 em 5 a partir de 5. E assim sucessivamente, até o último número da nossa lista. Os números não riscados ao final do método são números primos.

Vamos exemplificar. Suponha que listamos os números ímpares de 3 até 29.

3 5 7 9 11 13 15 17 19 21 23 25 27 29

A primeira etapa é pegar o primeiro número da lista, que é o 3, e riscar os números múltiplos de 3:

3 5 7 9 11 13 15 17 19 21 23 25 27 29

Pegamos o próximo número não riscado da lista, que é o 5, e riscamos os números múltiplos de 5:

3 5 7 9 11 13 15 17 19 21 23 25 27 29

O processo continua até acabarem os números não riscados. O resultado final do nosso exemplo é:

3 5 7 9 11 13 15 17 19 21 23 25 27 29

O que significa que os números 3, 5, 7, 11, 13, 17, 19, 23 e 29 são primos. Perceba que não realizamos contas nesse método, apenas riscamos os números múltiplos de um certo número em uma determinada etapa do método.

Agora, vamos ao que interessa: implementar o método!

Vamos começar definindo o vetor de números ímpares:


n = 29
v_num = seq(3,n,by=2)

> v_num
 [1]  3  5  7  9 11 13 15 17 19 21 23 25 27 29

Para sabermos quais números foram riscados e quais não foram, vamos definir um vetor do mesmo tamanho que v_num, cujas entradas sejam todas iguais a 1, ou seja, inicialmente não temos nenhum número riscado. O valor zero vai indicar que o número foi riscado.


tam_v = length(v_num)
v_i = rep(1,tam_v)
rbind(v_num,v_i)

> rbind(v_num,v_i)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
v_num    3    5    7    9   11   13   15   17   19    21    23    25    27    29
v_i      1    1    1    1    1    1    1    1    1     1     1     1     1     1

Desse modo, v_num[1] armazena o número 3 e v_i[1] armazena se o número 3 foi riscado (0) ou não (1).

Agora entramos na implementação do método:


Para cada elemento do vetor v_num
	pega o número primo atual
  	seja aux = indice atual + número primo atual
    	enquanto aux for menor que o tamanho do vetor v_num
    		altera para 0 o vetor v_i no indice aux
        	incrementa aux do valor do número primo atual
        fim enquanto
fim do para

O pseudo código acima pode ser escrito da seguinte forma no R:


for(i in 1:tam_v){
  p_atual = v_num[i];
  aux = i + p_atual;
  while(aux < tam_v){
    v_i[aux] = 0;
    aux = aux + p_atual;
  }
}

Ou seja, para cada elemento do vetor v_num, armazenamos o número primo atual em p_atual, armazenamos o próximo múltiplo desse número na variável aux e, enquanto não percorremos todo o vetor de números ímpares, vamos alterando para 0 o valor de v_i no índice correspondente aos múltiplos do primo atual. Ao final do processo, temos:


> rbind(v_num,v_i)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
v_num    3    5    7    9   11   13   15   17   19    21    23    25    27    29
v_i      1    1    1    0    1    1    0    1    1     0     1     0     0     1

Note que o resultado acima é exatamente o resultado do crivo de Eratóstenes: os números primos não estão riscados (1) e os números compostos estão riscados (0). Para selecionarmos apenas os número primos, fazemos:


> v_num[which(v_i==1)]
[1]  3  5  7 11 13 17 19 23 29

Como o código está solto, vamos colocá-lo dentro de uma função:


crivo_erat = function(n){
  v_num = seq(3,n,by=2)
  tam_v = length(v_num)
  v_i = rep(1,tam_v)
  
  for(i in 1:tam_v){
    p_atual = v_num[i];
    aux = i + p_atual;
    while(aux < tam_v){
      v_i[aux] = 0;
      aux = aux + p_atual;
    }
    
  }
  primos = v_num[which(v_i==1)];
  
  return(primos)
}

Desse modo, nossa função recebe um argumento: o número n escolhido, e retorna o resultado: um vetor contendo os números primos. Vamos testar:


> crivo_erat(100)
 [1]  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 99
 
> crivo_erat(1000)
  [1]   3   5   7  11  13  17  19  23  29  31  37  41  43  47  53  59  61  67  71  73  79  83
 [23]  89  97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197
 [45] 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331
 [67] 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461
 [89] 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607
[111] 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751
[133] 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907
[155] 911 919 929 937 941 947 953 967 971 977 983 991 997 999

Agora que temos o crivo de Eratóstenes implementado, podemos apontar dois pontos ineficientes da nossa implementação:

1) Alguns números são riscados da lista mais de uma vez.

2) Não é necessário percorrer todo o vetor de números ímpares para obtermos o resultado.

Vamos discutir esses dois pontos com o exemplo de uma lista de números de 3 a 49:

3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Primeiro ponto: na etapa do número primo 5, o número 15 já foi riscado na etapa do número primo 3. Note que 15 = 5*3, ou seja, o próximo número não riscado múltiplo de 5 do nosso vetor de números ímpares seria o número 5*5 = 25. Na etapa do número primo 7, o número 21 = 7*3 foi riscado por ser múltiplo de 3 e o número 35 = 7*5 foi riscado por ser múltiplo de 5. Ou seja, o próximo número não riscado múltiplo de 7 é 7*7 = 49. Dessa forma, dado um número primo atual p, podemos começar a riscar números a partir de p^2, evitando cortes desnecessários e tornando o algoritmo mais eficiente.

Segundo ponto: podemos encerrar o algoritmo quando o número primo atual for maior que a raiz quadrada do maior número da lista de ímpares n. Por quê? Lembrem das aulas passadas, quando otimizamos nosso algoritmo para encontrar números primos: um número composto terá um fator menor ou igual a sua raiz quadrada. Logo, se quisermos saber se 49 é número primo, basta verificarmos divisão até o número 7. Mas e os demais números do crivo? Esses números serão menores que o maior número da lista e, portanto, a verificação é ainda menor. Por exemplo, paramos a verificação do número 47 quando o primo atual é maior que 6,8. Isso significa que, depois do número 7, não cortaremos nada da lista de números ímpares.
Não te convenci? Vamos ao exemplo.
Após cortamos múltiplos de 3:

3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Após cortamos múltiplos de 5:

3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Após cortamos múltiplos de 7:

3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Agora, na etapa do número 11, o número 33 já foi cortado e não há mais múltiplos de 11 na nossa lista.

Na etapa do número 13, o número 39 já foi cortado e não há mais múltiplos de 13 na nossa lista.

Perceba que, se todos os múltiplos do número atual menores que seu quadrado já foram cortados, não há mais o que fazer nesse exemplo. No caso do número 11, começaríamos a cortar números a partir de 121, que não está na nossa lista. Dessa forma, podemos encerrar o algoritmo na etapa do número 7, que é a raiz quadrada de maior número da nossa lista.

Vamos implementar esses dois pontos no nosso algoritmo:


crivo_erat2 = function(n){
  v_num = seq(3,n,by=2)
  tam_v = length(v_num)
  v_i = rep(1,tam_v)
  ult_num = v_num[tam_v]
  
  for(i in 1:tam_v){
    if(v_num[i]>sqrt(ult_num)) break;
    p_atual = v_num[i];
    aux = (p_atual^2 - 1)/2;
    while(aux < tam_v){
      v_i[aux] = 0;
      aux = aux + p_atual;
    }
  }
  primos = v_num[which(v_i==1)];
  
  return(primos)
}

Vejamos as alterações que fizemos. O primeiro ponto está na seguinte linha de código


    aux = (p_atual^2 - 1)/2;

A linha de código acima diz o seguinte: o primeiro índice que pegaremos será o correspondente ao quadrado do número primo atual. Por exemplo, se estamos na etapa no número 7, começamos os cortes a partir do número 49, que tem como índice (49-1)/2 = 24 na nossa lista de ímpares.

O segundo ponto está implementado na seguinte linha de código


    if(v_num[i]>sqrt(ult_num)) break;

Simplesmente paramos o for quando o número primo atual é maior que a raiz quadrada do maior número da nossa lista.

Beleza, implementamos as melhorias, ou otimizações. Mas você deve estar se perguntando: vale a pena todo esse trabalho? É claro que a implementação não deu trabalho, são linhas de código simples, mas pensá-las é o que realmente dá trabalho e nem sempre é simples. Será que vale o esforço? Vamos analisar o custo de tempo das duas funções utilizando a função benchmark do pacote rbenchmark. Vamos aproveitar e comparar o tempo de execução com a função calcula_primos2() que implementamos nas aulas passadas.


> numero = 200;
> benchmark(crivo_erat(numero),crivo_erat2(numero),
+           calcula_primos2(numero),
+           replications = 1000,
+           columns = c("test", "replications", "elapsed","relative"))
                     test replications elapsed relative
3 calcula_primos2(numero)         1000    0.24      4.8
1      crivo_erat(numero)         1000    0.06      1.2
2     crivo_erat2(numero)         1000    0.05      1.0

> numero = 2000;
> benchmark(crivo_erat(numero),crivo_erat2(numero),
+           calcula_primos2(numero),
+           replications = 1000,
+           columns = c("test", "replications", "elapsed","relative"))
                     test replications elapsed relative
3 calcula_primos2(numero)         1000    3.06   16.105
1      crivo_erat(numero)         1000    0.56    2.947
2     crivo_erat2(numero)         1000    0.19    1.000

Percebam que sim! valeu a pena implementar as melhorias do crivo no nosso algoritmo. Vejam que mesmo para n pequeno, a função crivo_erat2 foi a de menor tempo. Para n = 2000, vemos uma grande diferença de tempo entre as duas funções do crivo, em que a função não otimizada leva quase 3 vezes o tempo de execução da função otimizada.

Agora, em comparação com a nossa função calcula_primos2 das aulas passadas (que já está otimizada), percebam que o crivo de Eratóstenes é bem mais eficiente. O que isso significa? Que apenas otimização de uma função não garante que ela seja a melhor, mas sim que ela é melhor que sua versão não otimizada. Podem existir métodos mais eficientes do qual você esteja usando, que valem mais a pena serem implementados do que simplesmente otimizar a função do método atual. Lembrem-se que, na prática, apesar de as funções (e métodos) retornarem o mesmo resultado, o melhor é escolher o mais eficiente!

Até a próxima aula!

Referência citada:

Coutinho, S.C. Números Inteiros e Criptografia RSA. 2005. Série de Computação e Matemática. IMPA.