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!