aula de hoje: redes de neurônios artificiais com o Keras

    1. introdução: o que é o Keras
    1. classificação dos dígitos do MNIST com redes totalmente conectadas
    1. classificação dos dígitos do MNIST com redes convolucionais
    1. regressão: redshifts fotométricos para o S-PLUS com redes totalmente conectadas

referência básica: Deep Learning with R, F. Chollet & J.J. Allaire

https://www.manning.com/books/deep-learning-with-r
segunda edição já disponível!!

https://github.com/jjallaire/deep-learning-with-r-notebooks/




1. introdução: o Keras

O Keras (https://keras.rstudio.com) é uma excelente plataforma para aplicações em deep learning (DL), já que permite definir e treinar quase qualquer tipo de modelo de DL. Ele provê interface para várias implementações de DL, como TensorFlow, Theano e CNTK.




2. classificação dos dígitos do MNIST com redes totalmente conectadas

Vamos ilustrar o Keras com um problema “clássico” de classificação de imagens: a classificação dos dígitos do MNIST

O MNIST é uma base de dados contendo imagens digitalizadas (28x28 pixels cada uma) dos dígitos 0 a 9, com 256 níveis de cinza.

Esta base de dados vem com o Keras. Vamos carregar os dados do MNIST e definir os conjuntos de treinamento e teste:

#tempo de processamento do script
t00 = Sys.time()

# reprodutibilidade
set.seed(1234)

# paralelização
library(doParallel)
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
cl <- makePSOCKcluster(5)
registerDoParallel(cl)

library(keras)
library(ggplot2)

# dados
mnist <- dataset_mnist()
## Loaded Tensorflow version 2.9.2
train_images <- mnist$train$x
train_labels <- mnist$train$y
test_images <- mnist$test$x
test_labels <- mnist$test$y

# Algo sobre os dados:
dim(train_images)
## [1] 60000    28    28
length(train_labels)
## [1] 60000
dim(test_images)
## [1] 10000    28    28
length(test_labels)
## [1] 10000
# os pixels das imagens tem valores entre 0 e 255
summary(train_images[,14,14])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00   29.00   96.97  232.00  255.00
# os targets são dígitos de 0 a 9
# exemplo:
test_labels[1:10]
##  [1] 7 2 1 0 4 1 4 9 5 9

Para visualizar as 4 primeiras imagens do conjunto de treinamento:

# visualização dos dígitos
par(mfcol=c(2,2))
par(mar=c(0, 0, 3, 0), xaxs='i', yaxs='i')
for (idx in 1:4) { 
    im <- train_images[idx,,]
    im <- t(apply(im, 2, rev)) 
    image(1:28, 1:28, im, col=gray((0:255)/255), 
          xaxt='n', main=paste(train_labels[idx]))
}


Vamos definir a arquitetura da rede: camadas, conexões e funções de ativação.

Vamos considerar uma rede densamente conectada com uma camada de entrada com 28x28=784 neurônios (parâmetro ‘input shape’), 512 numa camada oculta e 10 na camada de saída (uma saida para cada dígito): 784:512:10

A camada de saída vai dar a probabilidade de cada dígito, dado o input. Para isso vamos usar a softmax como função de ativação na saída: \[p_i = {e^{-y_i} \over \sum_k e^{-y_k}}\]

A definição de um modelo no Keras tem duas partes:

1.) as camadas e as funções de ativação- vamos considerar funções de ativação tipo ReLU na camada oculta e softmax na de saída:

# o "multilayer perceptron" (MLP) é um exemplo de rede densamente conectada 
# MLP com uma camada oculta com 512 unidades ReLU
modelo <- keras_model_sequential() %>%  
  layer_dense(units = 512, activation = "relu", input_shape = c(28 * 28)) %>% 
  layer_dense(units = 10, activation = "softmax")
  
# ativação 'softmax'  na saída: apropriada para classificação multi-classe

# se quiser incluir uma nova camada oculta densamente conectada com, por exemplo, 64 unidades e ativação tanh, inclua o seguinte comando entre a camada de entrada e a de saída:
# layer_dense(units = 64, activation = "tanh") %>%

2.) e os parâmetros de treinamento- otimização, função de custo e métrica de desempenho, adequados para um problema de classificação multi-classe:

lembrando
acurácia: fração de detecções corretas \[\frac{TP+TN}{TP+TN+FP+FN}\]

modelo %>% compile(
  optimizer = "rmsprop",                # variante turbinada da descida do gradiente
  loss = "categorical_crossentropy",    # classificação multi-classe
  metrics = c("accuracy")               # classificação
)

# um resumo do modelo:
summary(modelo)
## Model: "sequential"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  dense_1 (Dense)                    (None, 512)                     401920      
##  dense (Dense)                      (None, 10)                      5130        
## ================================================================================
## Total params: 407,050
## Trainable params: 407,050
## Non-trainable params: 0
## ________________________________________________________________________________

Vamos fazer um pré-processamento dos dados:

  1. transformando a imagem em um vetor
  1. colocando os valores dos pixels entre 0 e 1
# pré-processamento do conjunto de treino
dim(train_images)
## [1] 60000    28    28
# muda o formato do arquivo train_images para o keras entender
train_images <- array_reshape(train_images, c(60000, 28 * 28)) 

dim(train_images)
## [1] 60000   784
# agora cada imagem é representada por um vetor com 784 valores de pixels

# divide-se os valores dos pixels por 255 para ficarem entre 0 e 1
min(train_images)
## [1] 0
max(train_images)
## [1] 255
train_images <- train_images / 255  

# pré-processamento do conjunto de teste
test_images <- array_reshape(test_images, c(10000, 28 * 28)) 
test_images <- test_images / 255  

Este é um problema de classificação: as classes são variáveis categóricas

train_labels <- to_categorical(train_labels)
test_labels <- to_categorical(test_labels)

# as classes são 'hot-coded': todas 0, exceto uma (o target) que é 1:
head(train_labels)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    0    0    0    0    0    1    0    0    0     0
## [2,]    1    0    0    0    0    0    0    0    0     0
## [3,]    0    0    0    0    1    0    0    0    0     0
## [4,]    0    1    0    0    0    0    0    0    0     0
## [5,]    0    0    0    0    0    0    0    0    0     1
## [6,]    0    0    1    0    0    0    0    0    0     0

treinamento com 10 épocas:

épocas: número de vezes que a rede vê todo o conjunto de treunamento

batch-size: número de amostras que são propagadas por back-propagation antes que os pesos são atualizados

t0 = Sys.time()
historia <- modelo %>% fit(train_images, train_labels, epochs = 10, batch_size = 128,verbose=0,validation_data = list(test_images,test_labels))
plot(historia)

Sys.time()-t0
## Time difference of 25.67683 secs

note o baixo tempo de processamento: o Keras é muito eficiente!

desempenho do classificador no conjunto de teste:

metrics <- modelo %>% evaluate(test_images, test_labels, verbose = 0)
metrics
##       loss   accuracy 
## 0.06987681 0.98089999

esta rede tem um ótimo desempenho: a acurácia é de 98.4% no conjunto de teste




3. redes neurais convolucionais: CNNs ou convnets

As CNNs revolucionaram a análise de imagens e de certo modo viraram o padrão para este tipo de análise. Mas note que o algoritmo pode ser usado com muitos outros tipos de dados, desde que tenham um formato adequado.

O “coração” de uma CNN são as pilhas de camadas de convolução e pooling, layer_conv_2d() e layer_max_pooling_2d(), com uma ou duas camadas densamente conectadas na saída.

Uma camada de convolução passa um filtro de largura kernel_size sobre uma imagem ou feature map e transfere este sinal para um número grande de “filtros”, usando uma certa função de ativação, gerando um novo feature map.

A camada de max_pooling tem o papel de reduzir o tamanho dos mapas de filtros gerados pela camada de convolução.

Vamos ilustrar uma CNN com o problema da classificação dos dígitos do MNIST.

Uma imagem pode ser descrita como um tensor de forma (altura, largura, número de canais). As imagens do MNIST são (28,28,1).

# Vamos construir o modelo em duas etapas: a convolucional e depois a densamente conectada.

# parte convolucional do modelo
model <- keras_model_sequential() %>% 
  layer_conv_2d(filters = 32, kernel_size = c(3, 3), activation = "relu",
                input_shape = c(28, 28, 1)) %>% 
  layer_max_pooling_2d(pool_size = c(2, 2)) %>% 
  layer_conv_2d(filters = 64, kernel_size = c(3, 3), activation = "relu") %>% 
  layer_max_pooling_2d(pool_size = c(2, 2)) %>% 
  layer_conv_2d(filters = 64, kernel_size = c(3, 3), activation = "relu")

# resumo do modelo:
summary(model)
## Model: "sequential_1"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  conv2d_2 (Conv2D)                  (None, 26, 26, 32)              320         
##  max_pooling2d_1 (MaxPooling2D)     (None, 13, 13, 32)              0           
##  conv2d_1 (Conv2D)                  (None, 11, 11, 64)              18496       
##  max_pooling2d (MaxPooling2D)       (None, 5, 5, 64)                0           
##  conv2d (Conv2D)                    (None, 3, 3, 64)                36928       
## ================================================================================
## Total params: 55,744
## Trainable params: 55,744
## Non-trainable params: 0
## ________________________________________________________________________________
# notem a estrutura de funil

A saida de cada camada layer_conv_2d() e layer_max_pooling_2d() é um tensor tri-dimensional de forma (altura,largura,canais). A largura e a altura tendem a diminuir com a profundidade na rede.

O número de canais é controlado pelo primeiro argumento do layer_conv_2d() (32 ou 64)

A próxima etapa é conectar a saída desse modelo, um tensor de forma (3,3,64) a uma rede densamente conectada para classificação.

A entrada dessas redes são vetores, de modo que é necessário transformar o tensor 3D em um vetor 1D; o nome dessa operação no Keras é “flatten¨.

Vamos considerar uma rede densa de apenas uma camada:

model <- model %>% 
# flatten
  layer_flatten() %>% 
# camada densamente conectada  
  layer_dense(units = 64, activation = "relu") %>%   
# saída - classificação 
  layer_dense(units = 10, activation = "softmax")     

# resumo do modelo:
summary(model) 
## Model: "sequential_1"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  conv2d_2 (Conv2D)                  (None, 26, 26, 32)              320         
##  max_pooling2d_1 (MaxPooling2D)     (None, 13, 13, 32)              0           
##  conv2d_1 (Conv2D)                  (None, 11, 11, 64)              18496       
##  max_pooling2d (MaxPooling2D)       (None, 5, 5, 64)                0           
##  conv2d (Conv2D)                    (None, 3, 3, 64)                36928       
##  flatten (Flatten)                  (None, 576)                     0           
##  dense_3 (Dense)                    (None, 64)                      36928       
##  dense_2 (Dense)                    (None, 10)                      650         
## ================================================================================
## Total params: 93,322
## Trainable params: 93,322
## Non-trainable params: 0
## ________________________________________________________________________________
# o tensor (3, 3, 64) é "achatado" em um vetor de dimensão 3*3*64 = 576 antes de passar pelas camadas densas.

Vamos carregar os dados do MNIST, na forma de conjuntos de treino e teste e num formato adequado para este tipo de rede:

# para monitorar o tempo de processamento deste módulo
t0 = Sys.time()  

mnist <- dataset_mnist()
c(c(train_images, train_labels), c(test_images, test_labels)) %<-% mnist
dim(train_images)
## [1] 60000    28    28
summary(train_images[,14,14])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00   29.00   96.97  232.00  255.00
# vamos transformar o tensor num vetor
train_images <- array_reshape(train_images, c(60000, 28, 28, 1))

# vamos colocar os valores dos pixels entre 0 e 1:
train_images <- train_images / 255

test_images <- array_reshape(test_images, c(10000, 28, 28, 1))
test_images <- test_images / 255

# em classificação os labels são variáveis categóricas:
train_labels <- to_categorical(train_labels)
test_labels <- to_categorical(test_labels)

# hiper-parâmetros do modelo
model %>% compile(
  optimizer = "rmsprop",
  loss = "categorical_crossentropy",
  metrics = c("accuracy")
)

# ajuste            
model %>% fit(
  train_images, train_labels, 
  epochs = 5, batch_size=64,
  verbose=0
)

# avaliação dos resultados no conjunto de teste
results <- model %>% evaluate(test_images, test_labels)

results
##       loss   accuracy 
## 0.03173341 0.99070001
Sys.time() - t0
## Time difference of 1.993904 mins

compare a acurácia do conjunto de teste com o valor obtido com a rede densa anterior.



4. regressão: redshift fotométrico com a fotometria do S-PLUS com redes totalmente conectadas

Vamos usar o mesmo conjunto de dados da aula anterior. Vamos carregar a biblioteca caret para usar suas rotinas de partição dos dados em treinamento e teste.

library(caret)
## Loading required package: lattice
dados = as.data.frame(read.table('splus-mag-z.dat', sep = "", header=TRUE))

# vamos examinar os dados:
dim(dados)
## [1] 55803    13
# cada linha 
head(dados)
##   uJAVA_petro F378_petro F395_petro F410_petro F430_petro g_petro F515_petro
## 1       18.71      18.69      18.50      18.12      17.78   17.43      17.09
## 2       18.26      18.12      18.28      17.32      16.95   16.67      16.24
## 3       19.38      19.45      18.76      18.60      18.50   18.07      17.69
## 4       20.30      19.73      20.01      19.18      18.91   18.36      18.00
## 5       19.49      18.98      18.64      18.64      18.02   17.14      16.61
## 6       19.86      19.38      18.93      18.75      18.15   17.44      16.94
##   r_petro F660_petro i_petro F861_petro z_petro z_SDSS
## 1   16.89      16.83   16.56      16.54   16.45  0.111
## 2   15.97      15.90   15.61      15.52   15.41  0.082
## 3   17.55      17.47   17.25      17.31   17.18  0.086
## 4   17.75      17.68   17.51      17.46   17.34  0.111
## 5   16.06      15.93   15.65      15.54   15.38  0.162
## 6   16.58      16.48   16.24      16.11   16.00  0.082
# para o tempo de processamento não ficar muito grande, vamos considerar uma amostra com 10000 galáxias
ngal = 10000
dados = dados[1:ngal,]
z_SDSS = dados[,13]

# número das linhas dos objetos selecionados para o conjunto de treinamento
numlinhastreinamento =  createDataPartition(z_SDSS, p=0.75, list=FALSE)

# conjunto de treinamento:
treinamento =  dados[numlinhastreinamento,]

# conjunto de teste
teste = dados[-numlinhastreinamento,]

train_data = treinamento[,1:12]
train_targets = treinamento[,13]

test_data = teste[,1:12]
test_targets = teste[,13]

pré-processamento dos dados:

vamos padronizar os dados de entrada, as magnitudes (subtração da média e divisão pelo desvio padrão):

mean <- apply(train_data, 2, mean)
std <- apply(train_data, 2, sd)
train_data <- scale(train_data, center = mean, scale = std)
test_data <- scale(test_data, center = mean, scale = std)

arquitetura da rede- duas camadas densas (com função de ativação ReLU) e uma unidade de saída linear (apropriada para regressão):

# Como vamos usar o mesmo modelo mais de uma vez, vamos usar uma função para construí-lo:
build_model <- function() {
  model <- keras_model_sequential() %>% 
    layer_dense(units = 32, activation = "relu", 
                input_shape = dim(train_data)[[2]]) %>% 
    layer_dense(units = 32, activation = "relu") %>% 
    layer_dense(units = 1) #no activation = linear layer
    
  model %>% compile(
    optimizer = "rmsprop", 
    loss = "mse",             #Mean Squared Error
    metrics = c("mae")        #Mean Absolute Error
  )
}

Vamos treinar inicialmente o modelo com 1000 épocas:

t0 = Sys.time()

  model <- build_model()
  historia <- model %>% fit(train_data, train_targets,
                epochs = 1000, batch_size = 16, verbose = 0,validation_data = list(test_data,test_targets))
                
 summary(model)               
## Model: "sequential_2"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  dense_6 (Dense)                    (None, 32)                      416         
##  dense_5 (Dense)                    (None, 32)                      1056        
##  dense_4 (Dense)                    (None, 1)                       33          
## ================================================================================
## Total params: 1,505
## Trainable params: 1,505
## Non-trainable params: 0
## ________________________________________________________________________________
 plot(historia)     

 Sys.time()-t0
## Time difference of 13.14116 mins

Com 1000 épocas temos um claro overfitting

Vamos treinar o modelo com 100 épocas e verificar sua convergência e se houve overfitting:

historia <- model %>% fit(train_data, train_targets, epochs = 100, batch_size = 128,verbose=0,validation_data = list(test_data,test_targets))              
 
 plot(historia)

results <- model %>% evaluate(test_data, test_targets, verbose = 0)
results
##        loss         mae 
## 0.003480751 0.039228413
# predições

z_pred = model %>% predict(test_data)
z_pred = z_pred[,1]

# estatística da performance do algoritmo
# sigma equivalente da gaussiana
relz = (z_pred-test_targets)
sig_G = 0.7413*(quantile(relz,0.75,names = FALSE) - quantile(relz,0.25,names = FALSE))
sig_G
## [1] 0.03697391
# na aula anterior obtivemos 0.035 e 0.040 com rf e knn, respectivamente

#visualização do resultado
my_data = as.data.frame(cbind(predicted = z_pred,
                            observed = test_targets))
# Plot predictions vs test data
ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = "darkred", alpha = 0.5) + 
      geom_smooth(method=lm)+ ggtitle('GLM ') +
      ggtitle("DL: prediction vs test data") +
      xlab("prediction ") +
      ylab("observed data") +
      theme(plot.title = element_text(color="darkgreen",size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
            axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))
## `geom_smooth()` using formula 'y ~ x'

Sys.time()-t0
## Time difference of 13.354 mins

tempo de processamento do script:

Sys.time()-t00
## Time difference of 15.94609 mins
# sem paralelizaçã0: Time difference of 11.41628 mins

Exercício

Examine como a escolha da arquitetura do modelo afeta a acurácia do conjunto de teste no caso da classificação dos dígitos do MNIST. Considere o exemplo de nossa rede densamente conectada. No exemplo, ela tem 512 neurônios em uma camada oculta.

    1. inicialize a semente aleatória com seu número USP
    1. varie o número de neurônios em cada camada: use 32 e 64
    1. varie o número de camadas ocultas (com 64 neurônios em cada uma): use de 1 a 3
    1. varie a ativação no modelo com uma camada com 512 neurônios: compare ReLU com tanh

É bom variar um pouco o número de épocas para encontrar um valor que evita um grande overfitting ou underfitting (não precisa otimizar)

Escreva uma avaliação sobre o desempenho da acurácia do modelo em função desses parâmetros.