# ======================================================================================================== # # STEPHANE TUFFERY # BIG DATA, MACHINE LEARNING ET APPRENTISSAGE PROFOND # EDITIONS TECHNIP # ======================================================================================================== # # --------------------------------------------------------------------------------------------------------- # Chapitre 1 : Le Big Data # --------------------------------------------------------------------------------------------------------- # Section 1.3.1 : Prédictions avec Google Trends ("à la Flu Trends") # --------------------------------------------------------------------------------------------------------- # CAC 40 avec quantmod install.packages('quantmod') library(quantmod) getSymbols("^FCHI", src="yahoo") head(FCHI) plot(FCHI[,c("FCHI.Adjusted")], type="l") chartSeries(FCHI["2007::"]) cac40_cours <- as.data.frame(FCHI) # suppression du 1/5/2014 dont les valeurs sont manquantes cac40_cours <- cac40_cours[which(!is.na(cac40_cours$FCHI.Volume)),] cac40_cours <- cbind("date" = rownames(cac40_cours), cac40_cours) rownames(cac40_cours) <- NULL cac40_cours$mois <- substr(cac40_cours$date,1,7) head(cac40_cours) # calcul du cours de clôture ajusté mensuel moyen cac40_mensuel <- aggregate(cac40_cours[,"FCHI.Adjusted"], list(date=cac40_cours$mois), mean) # mise de la date sous forme 01-mois-année de la classe habituelle POSIXct cac40_mensuel$date <- as.POSIXct(paste(cac40_mensuel$date, "01", sep="-")) colnames(cac40_mensuel) <- c("date", "cours") head(cac40_mensuel) # Google Trends install.packages('gtrendsR') library(gtrendsR) cac40_trend <- gtrends("CAC 40", geo="FR", time="2007-01-01 2015-12-31") #cac40_trend <- gtrends("CAC 40", geo="FR", time = "today+5-y") cac40_trend <- gtrends("CAC 40", geo="FR", time = "all") # avec une résolution “1h”, “4h”, “1d”, “7d”, Google Trends renvoie les données de la dernière période d'1h, 4h, 1 jour et 7 jours # et les paramètres start_date et end_date sont inutiles # une résolution non indiquée est mensuelle au-delà de 5 ans, hebdomadaire pour toute période entre 3 mois et 5 ans, # et quotidienne pour toute période entre 1 et 3 mois head(cac40_trend$interest_over_time) plot(cac40_trend) # comparaison des cours du CAC40 et des tendances de Google cac <- merge(cac40_mensuel, cac40_trend$interest_over_time[,1:2]) head(cac) tail(cac) # corrélations cor(cac$cours, cac$hits) cor(cac$cours, cac$hits, method="spearman") plot(cac$cours, cac$hits) abline(lm(cac$hits ~ cac$cours)) # séries temporelles cac.ts <- ts(data = cac, start = c(2007,1), frequency = 12) cac.ts plot(cac.ts[,-1]) plot(scale(cac.ts[,-1]), plot.type = "single", lty = 1:2) # variables centrées-réduites pour mise sur une même échelle # cross-corrélation (de Pearson) cross <- ccf(cac.ts[,2],cac.ts[,3], type = c("correlation"), plot = TRUE) # on voit que la corrélation linéaire est maximale quand la 2e série (Google trends) précède la 1ère série (cours de Bourse) de 1 mois cross$lag[which.max(abs(cross$acf))]*12 # nb de mois de décalage maximisant la cross-corrélation en valeur absolue cross$acf[which.max(abs(cross$acf))] # cross-corrélation maximale en valeur absolue # on voit que la corrélation des rangs est maximale quand la 2e série (Google trends) précède la 1ère série (cours de Bourse) de 1 mois cross <- ccf(rank(cac.ts[,2]),rank(cac.ts[,3]), type = c("correlation"), plot = TRUE) cross$lag[which.max(abs(cross$acf))] # nb de mois de décalage maximisant la cross-corrélation en valeur absolue cross$acf[which.max(abs(cross$acf))] # cross-corrélation maximale en valeur absolue # corrélation par année library(lubridate) # pour sa fonction "year" cac$an <- year(cac$date) sapply(2007:2015, function(x) {cor(cac[which(cac$an==x),]$cours, cac[which(cac$an==x),]$hits, method="spearman")}) sapply(2007:2015, function(x) {scac <- window(cac.ts, start=c(x,1), end=c(x,12)) ; cor(scac[,2],scac[,3], method="spearman")}) # ajout de nouvelles recherches cac40_trend_tout <- gtrends(c("crise","dette","chômage","inflation","croissance"), geo="FR", time="2007-01-01 2015-12-31") head(cac40_trend_tout$interest_over_time) plot(cac40_trend_tout) CAC <- merge(cac40_mensuel, cac40_trend_tout$interest_over_time[,1:3]) sapply(c("crise","dette","chômage","inflation","croissance"), function(x) {cor(CAC[which(CAC$keyword==x),]$cours, CAC[which(CAC$keyword==x),]$hits, method="spearman")}) # --------------------------------------------------------------------------------------------------------- # Chapitre 2 : Les outils informatiques pour le Big Data # --------------------------------------------------------------------------------------------------------- # Section 2.6 : Web Scraping # --------------------------------------------------------------------------------------------------------- library(rvest) url <- "https://en.wikipedia.org/wiki/List_of_countries_and_dependencies_by_population" pop <- read_html(url) %>% html_table(fill=TRUE) %>% .[[1]] # --------------------------------------------------------------------------------------------------------- # Chapitre 5 : L'apprentissage profond (Deep Learning) # --------------------------------------------------------------------------------------------------------- # Section 5.3 : Perceptron # --------------------------------------------------------------------------------------------------------- sigmoid <- function(x) {1 / (1 + exp(-x))} feedforward <- function(x, w1, w2) { z1 <- cbind(1, x) %*% w1 h <- sigmoid(z1) z2 <- cbind(1, h) %*% w2 list(output = sigmoid(z2), h = h) } backpropagate <- function(x, y, y_hat, w1, w2, h, learn_rate) { dw2 <- t(cbind(1, h)) %*% (y_hat - y) %*% y %*% (1-y) dh <- (y_hat - y) %*% t(w2[-1,]) dw1 <- t(cbind(1, x)) %*% (h * (1 - h) * dh) w1 <- w1 - learn_rate * dw1 w2 <- w2 - learn_rate * dw2 list(w1 = w1, w2 = w2) } perceptron <- function(x, y, hidden = 5, learn_rate = 1e-2, iterations = 1e4) { d <- ncol(x) + 1 set.seed(235) w1 <- matrix(rnorm(d * hidden), d, hidden) w2 <- as.matrix(rnorm(hidden + 1)) for (i in 1:iterations) { ff <- feedforward(x, w1, w2) bp <- backpropagate(x, y, y_hat = ff$output, w1, w2, h = ff$h, learn_rate = learn_rate) w1 <- bp$w1; w2 <- bp$w2 } list(output = ff$output, w1 = w1, w2 = w2) } # perceptron sur jeu de données spam library(kernlab) data(spam) summary(spam) x <- data.matrix(spam[,-58]) set.seed(235) x <- scale(data.matrix(spam[sample(nrow(spam)),-58])) # normalisation + shuffle y <- spam$type == "spam" system.time(mynnet <- perceptron(x, y, hidden = 100, iterations = 1e3)) summary(mynnet$output) # évaluation du modèle library(pROC) auc(y, as.numeric(mynnet$output), quiet=TRUE) plot.roc(y, as.numeric(mynnet$output)) (tab <- table(y, (mynnet$output>0.5))) sum(diag(tab)/nrow(spam)) cbind(prop.table(tab,1), addmargins(tab,2)) # Section 5.14 : Algorithmes adaptatifs # --------------------------------------------------------------------------------------------------------- # installation mxnet : https://mxnet.incubator.apache.org/install/index.html # installation mxnet : https://mxnet.apache.org/get_started/ cran <- getOption("repos") cran["dmlc"] <- "https://apache-mxnet.s3-accelerate.dualstack.amazonaws.com/R/CRAN/" options(repos = cran) install.packages("mxnet") # version de mxnet pour R >= 3.6.0 install.packages("https://s3.ca-central-1.amazonaws.com/jeremiedb/share/mxnet/CPU/3.6/mxnet.zip", repos = NULL) # données directement disponibles au format csv : http://www.pjreddie.com/projects/mnist-in-csv/ train <- read.csv("D:/Data/MNIST/mnist_train.csv", header=TRUE) test <- read.csv("D:/Data/MNIST/mnist_test.csv", header=TRUE) library(mxnet) train.x <- train[,-1] train.y <- train[,1] train.x <- t(train.x/255) test.x <- test[,-1] test.x <- t(test.x/255) # 78 # transformation des matrices en arrays (nécessaire pour l'opérateur de convolution de mxnet) : train.array <- train.x dim(train.array) <- c(28, 28, 1, ncol(train.x)) test.array <- test.x dim(test.array) <- c(28, 28, 1, ncol(test.x)) # 2 couches de convolution avec pooling recouvrant et dropout # entrée data <- mx.symbol.Variable('data') # 1ère couche de convolution conv1 <- mx.symbol.Convolution(data=data, kernel=c(5,5), num_filter=20, workspace=2048) tanh1 <- mx.symbol.Activation(data=conv1, act_type="relu") pool1 <- mx.symbol.Pooling(data=tanh1, pool_type="max", kernel=c(3,3), stride=c(2,2)) # 2e couche de convolution conv2 <- mx.symbol.Convolution(data=pool1, kernel=c(5,5), num_filter=60, workspace=4096) tanh2 <- mx.symbol.Activation(data=conv2, act_type="relu") pool2 <- mx.symbol.Pooling(data=tanh2, pool_type="max", kernel=c(3,3), stride=c(2,2)) # couche complètement connectée flatten <- mx.symbol.Flatten(data=pool2) drp <- mx.symbol.Dropout(data=flatten, p=0.5) fc1 <- mx.symbol.FullyConnected(data=drp, num_hidden=500) tanh3 <- mx.symbol.Activation(data=fc1, act_type="relu") # couche de sortie fc2 <- mx.symbol.FullyConnected(data=tanh3, num_hidden=10) lenet <- mx.symbol.SoftmaxOutput(data=fc2) device.cpu <- lapply(0:7, function(i) { mx.cpu(i) }) system.time(model <- mx.model.FeedForward.create(symbol=lenet, X=train.array, y=train.y, ctx=device.cpu, num.round=10, array.batch.size=128, optimizer="sgd", learning.rate=0.1, momentum=0.65, eval.data=list(data=test.array, label=test$label), wd=0.00001, eval.metric = mx.metric.accuracy)) system.time(model <- mx.model.FeedForward.create(symbol=lenet, X=train.array, y=train.y, ctx=device.cpu, num.round=10, array.batch.size=128, optimizer="adagrad", eval.data=list(data=test.array, label=test$label), wd=0.00001, eval.metric = mx.metric.accuracy)) system.time(model <- mx.model.FeedForward.create(symbol=lenet, X=train.array, y=train.y, ctx=device.cpu, num.round=10, array.batch.size=128, optimizer="adadelta", eval.data=list(data=test.array, label=test$label), wd=0.00001, eval.metric = mx.metric.accuracy)) system.time(model <- mx.model.FeedForward.create(symbol=lenet, X=train.array, y=train.y, ctx=device.cpu, num.round=50, array.batch.size=128, optimizer="adadelta", eval.data=list(data=test.array, label=test$label), wd=0.00001, eval.metric = mx.metric.accuracy)) system.time(model <- mx.model.FeedForward.create(symbol=lenet, X=train.array, y=train.y, ctx=device.cpu, num.round=50, array.batch.size=128, optimizer="rmsprop", eval.data=list(data=test.array, label=test$label), wd=0.001, eval.metric = mx.metric.accuracy)) # Section 5.18 : Auto-encodeurs # --------------------------------------------------------------------------------------------------------- # installation de Keras à partir du CRAN install.packages('keras') # installation de Keras à partir de GitHub install.packages("devtools") devtools::install_github("rstudio/keras") library(keras) # ligne suivante à ne lancer que la première fois install_keras() # installe aussi tensorflow dont keras est une interface # préalablement : installer Anaconda pour Python 3.6 (https://www.continuum.io/downloads#windows) # en décembre 2019 : # - Keras n'est pas compatible avec Python 3.7 # - le package keras de R nécessite la version 1.12 du package reticulate et non la version plus récente 1.13 ou 1.14 # (c'est vrai que l'on soit sous R 3.5.3 ou R 3.6.1 : on peut charger le package keras et définir un modèle de réseau de neurones, # mais la fonction "fit" d'apprentissage se plante) # la récupération de la version 1.12 peut se faire comme suit : require(devtools) install_version("reticulate", version = "1.12", repos = "http://cran.us.r-project.org") library(keras) # chargement du jeu de données c(c(x_train, y_train), c(x_test, y_test)) %<-% dataset_mnist() # autre mode de chargement mnist <- dataset_mnist() x_train <- mnist$train$x y_train <- mnist$train$y x_test <- mnist$test$x y_test <- mnist$test$y # dimensions des images en entrée img_rows <- 28 img_cols <- 28 num_classes <- 10 # conversion d'un tenseur 3D (60000 x 28 x 28) en un tenseur 2D (60000 x 784) x_train <- array(x_train, dim = c(nrow(x_train), img_rows*img_cols)) x_test <- array(x_test, dim = c(nrow(x_test), img_rows*img_cols)) model <- keras_model_sequential() model %>% layer_dense(units = 400, activation = 'tanh', input_shape = ncol(x_train)) %>% layer_dense(units = 200, activation = 'tanh') %>% layer_dense(units = 100, activation = 'tanh') %>% layer_dense(units = 2, activation = 'tanh') %>% layer_dense(units = 100, activation = 'tanh') %>% layer_dense(units = 200, activation = 'tanh') %>% layer_dense(units = 400, activation = 'tanh') %>% layer_dense(units = ncol(x_train)) model # apprentissage de l'autoencodeur sur l'échantillon d'apprentissage model %>% compile(loss = "mean_squared_error", optimizer = "adamax") system.time(history <- fit(model, x_train, x_train, epochs = 100, batch_size = 64, validation_data = list(x_test, x_test), verbose=1)) plot(history) layer_name <- 'dense_3' intermediate_layer_model <- keras_model(inputs = model$input, outputs = get_layer(model, layer_name)$output) intermediate_output <- predict(intermediate_layer_model, x_train) head(intermediate_output) summary(intermediate_output) library(ggplot2) qplot(intermediate_output[,1], intermediate_output[,2], colour = as.factor(y_train)) # autoencodage avec une forêt aléatoire non-supervisée train <- read.csv("D:/Data/MNIST/mnist_train.csv", header=TRUE) test <- read.csv("D:/Data/MNIST/mnist_test.csv", header=TRUE) library(randomForest) set.seed(235) system.time(rf <- randomForest(train[,-1], ntree=500, mtry=30, nodesize=5, proximity=TRUE, keep.forest=F)) system.time(rf <- randomForest(test[,-1], ntree=500, mtry=30, nodesize=5, proximity=TRUE, keep.forest=F)) cible <- as.factor(test[,1]) MDSplot(rf, cible, palette=rainbow(10), pch=as.numeric(cible)) legend("topleft", levels(cible), col=rainbow(10), cex=0.8, fill=rainbow(10)) # autoencodage avec une forêt aléatoire supervisée library(ranger) train$label <- as.factor(train[,1]) system.time(rf <- ranger(label ~ ., data = train, num.trees=300, mtry=30, write.forest = TRUE, min.node.size=5, seed=235)) # extraction de la matrice de proximité install.packages('edarf') library(edarf) prox <- extract_proximity(rf, test) pca <- princomp(prox, cor = TRUE) plot_prox(pca, color = as.factor(test$label), color_label = "Label", size = 2) # analyse en composantes principales ACP <- princomp(train[,-1]) library(ggplot2) qplot(ACP$scores[,1], ACP$scores[,2], colour=as.factor(train[,1])) # Section 5.19.1 : YOLO # --------------------------------------------------------------------------------------------------------- #install.packages("devtools") si devtools n'est pas installé #devtools::install_github("bnosac/image", subdir = "image.darknet", build_vignettes = TRUE) library(image.darknet) param <- system.file(package="image.darknet", "include", "darknet", "cfg", "tiny-yolo-voc.cfg") poids <- system.file(package="image.darknet", "models", "tiny-yolo-voc.weights") lname = system.file(package = "image.darknet", "include", "darknet", "data", "voc.names") param <- system.file(package="image.darknet", "include", "darknet", "cfg", "yolov2.cfg") poids <- system.file(package="image.darknet", "models", "yolov2.weights") lname = system.file(package = "image.darknet", "include", "darknet", "data", "coco.names") yolo <- image_darknet_model(type = 'detect', model = param, weights = poids, labels = lname) f <- system.file("include", "darknet", "data", "dog.jpg", package="image.darknet") x <- image_darknet_detect(file = f, object = yolo) x <- image_darknet_detect(file = "D:/Data/Images/chatchien-1080x506.jpg", object = yolo, threshold = 0.05) # Section 5.19.5 : Analyse du style d’une œuvre picturale # --------------------------------------------------------------------------------------------------------- # Référence : François Chollet and J.J. Allaire (2018). "Deep Learning with R", section 8.3.3. # --------------------------------------------------------------------------------------------------------- # Chapitre 6 : Apprentissage profond en pratique # --------------------------------------------------------------------------------------------------------- # Installation du package mxnet : voir section 5.14 # Section 6.2.1 : Généralités sur MXNet # --------------------------------------------------------------------------------------------------------- library(mxnet) train <- read.csv("D:/Data/MNIST/mnist_train.csv", header=TRUE) test <- read.csv("D:/Data/MNIST/mnist_test.csv", header=TRUE) # préparation des données pour mxnet train.x <- train[,-1] train.y <- train[,1] train.x <- t(train.x/255) test.x <- test[,-1] test.x <- t(test.x/255) # 784 lignes et 10000 colonnes, valeurs entre 0 et 1 # 1er réseau : perceptron à 1 couche cachée mx.set.seed(0) model <- mx.mlp(train.x, train.y, hidden_node=100, out_node=10, out_activation="softmax", num.round=20, array.batch.size=128, learning.rate=0.07, momentum=0.9, eval.metric = mx.metric.accuracy) model <- mx.mlp(train.x, train.y, hidden_node=100, out_node=10, out_activation="softmax", num.round=20, array.batch.size=128, optimizer="adadelta", eval.metric = mx.metric.accuracy) # 2e réseau : perceptron à 2 couches cachées mx.set.seed(0) model <- mx.mlp(train.x, train.y, hidden_node=c(128,64), out_node=10, out_activation="softmax", num.round=10, array.batch.size=128, optimizer="adadelta", eval.metric = mx.metric.accuracy) # avec dropout model <- mx.mlp(train.x, train.y, hidden_node=c(128,64), out_node=10, out_activation="softmax", num.round=10, dropout=0.5, array.batch.size=128, optimizer="adadelta", eval.metric = mx.metric.accuracy) # prédiction preds <- predict(model, test.x) dim(preds) # 1 ligne par valeur à prédire - 1 colonne par observation # après transposition, on a une ligne par observation, et une colonne par valeur de 0 à 9 (la colonne 1 correspond à 0...) pred.label <- max.col(t(preds)) - 1 (nrow(test)-sum(test$label == pred.label))/nrow(test) # taux d'erreur table(test$label,pred.label) # matrice de confusion # autre façon de spécifier un perceptron à 1 couche cachée data <- mx.symbol.Variable("data") fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=100) act1 <- mx.symbol.Activation(fc1, name="relu1", act_type="relu") fc2 <- mx.symbol.FullyConnected(act1, name="fc3", num_hidden=10) out <- mx.symbol.SoftmaxOutput(fc2, name="sm") # autre façon de spécifier un perceptron multi-couches (2 couches cachées) data <- mx.symbol.Variable("data") fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=128) act1 <- mx.symbol.Activation(fc1, name="relu1", act_type="relu") fc2 <- mx.symbol.FullyConnected(act1, name="fc2", num_hidden=64) act2 <- mx.symbol.Activation(fc2, name="relu2", act_type="relu") fc3 <- mx.symbol.FullyConnected(act2, name="fc3", num_hidden=10) out <- mx.symbol.SoftmaxOutput(fc3, name="sm") # apprentissage du réseau mx.set.seed(0) ptm <- proc.time() model <- mx.model.FeedForward.create(softmax, X=train.x, y=train.y, ctx=mx.cpu(), num.round=20, array.batch.size=128, learning.rate=0.07, momentum=0.9, eval.metric=mx.metric.accuracy, initializer=mx.init.uniform(0.07), epoch.end.callback=mx.callback.log.train.metric(100)) proc.time() - ptm # Section 6.2.2 : Création d’un réseau convolutif avec MXNet # --------------------------------------------------------------------------------------------------------- # affichage des 100 premiers chiffres de MNIST old <- par(no.readonly = TRUE) par(mfrow = c(10,10), mai = c(0,0,0,0)) for (i in (1:100)) { im <- matrix(data=as.numeric((train[i,-1]>0)), nrow=28, ncol=28) rim <- matrix(NA,28,28) for (j in 1:28) {rim[,j] <- im[,29-j]} image(1:28, 1:28, rim, col=gray((0:255)/255), xlab=" ", ylab=" ", axes=F) text(3, 0, train.y[i], cex = 2, col = gray(0.75), font=3, pos = c(3,4)) title(main = train[i,1]) } par(old) # transformation des matrices en arrays (nécessaire pour l'opérateur de convolution de mxnet) train.array <- train.x dim(train.array) <- c(28, 28, 1, ncol(train.x)) test.array <- test.x dim(test.array) <- c(28, 28, 1, ncol(test.x)) # réseau à 2 couches de convolution # entrée data <- mx.symbol.Variable('data') # première couche de convolution conv1 <- mx.symbol.Convolution(data=data, kernel=c(5,5), num_filter=20) tanh1 <- mx.symbol.Activation(data=conv1, act_type="tanh") pool1 <- mx.symbol.Pooling(data=tanh1, pool_type="max", kernel=c(2,2), stride=c(2,2)) # deuxième couche de convolution conv2 <- mx.symbol.Convolution(data=pool1, kernel=c(5,5), num_filter=50) tanh2 <- mx.symbol.Activation(data=conv2, act_type="tanh") pool2 <- mx.symbol.Pooling(data=tanh2, pool_type="max", kernel=c(2,2), stride=c(2,2)) # première couche complètement connectée flatten <- mx.symbol.Flatten(data=pool2) fc1 <- mx.symbol.FullyConnected(data=flatten, num_hidden=500) tanh3 <- mx.symbol.Activation(data=fc1, act_type="tanh") # deuxième couche complètement connectée fc2 <- mx.symbol.FullyConnected(data=tanh3, num_hidden=10) # sortie lenet <- mx.symbol.SoftmaxOutput(data=fc2) # --------------------------------------------------------------------------------------- # 2 couches de convolution avec normalisation batch # input library(dplyr) lenet <- mx.symbol.Variable('data') %>% # première couche de convolution mx.symbol.Convolution(kernel=c(5,5), num_filter=20) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(2,2), stride=c(2,2)) %>% #mx.symbol.Dropout(p=0.25) %>% # deuxième couche de convolution mx.symbol.Convolution(kernel=c(5,5), num_filter=50) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(2,2), stride=c(2,2)) %>% #mx.symbol.Dropout(p=0.25) %>% # première couche complètement connectée mx.symbol.Flatten() %>% mx.symbol.FullyConnected(num_hidden=500) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # deuxième couche complètement connectée mx.symbol.FullyConnected(num_hidden=10) %>% # sortie mx.symbol.SoftmaxOutput() # --------------------------------------------------------------------------------------- # 2 couches de convolution avec normalisation batch et dropout # entrée lenet <- mx.symbol.Variable('data') %>% # première couche de convolution mx.symbol.Convolution(kernel=c(5,5), num_filter=20) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(3,3), stride=c(2,2)) %>% # deuxième couche de convolution mx.symbol.Convolution(kernel=c(5,5), num_filter=60) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(3,3), stride=c(2,2)) %>% # première couche complètement connectée mx.symbol.Flatten() %>% mx.symbol.Dropout(p=0.1) %>% mx.symbol.FullyConnected(num_hidden=600) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # deuxième couche complètement connectée mx.symbol.FullyConnected(num_hidden=10) %>% # sortie mx.symbol.SoftmaxOutput() # apprentissage sur toute la CPU disponible device <- mx.cpu() # apprentissage logger <- mx.metric.logger$new() mx.set.seed(235) temps <- proc.time() model <- mx.model.FeedForward.create(lenet, X=train.array, y=train.y, ctx=device, num.round=30, array.batch.size=128, optimizer="adadelta", wd=0.00001, eval.data=list(data=test.array, label=test$label), eval.metric=mx.metric.accuracy, epoch.end.callback=mx.callback.log.train.metric(128,logger)) proc.time() - temps # courbe de l'historique du taux de bon classement plot(logger$eval, type="l", ylim=c(0.9,1), ylab="taux bon classement", xlab="epoch") abline(h=logger$eval[which.max(logger$eval)], lty=2) lines(logger$train, lty=3) # application du modèle preds <- predict(model, test.array) pred.label <- max.col(t(preds)) - 1 (nrow(test)-sum(test$label == pred.label))/nrow(test) # taux d'erreur table(test$label,pred.label) # représentation du réseau graph.viz(model$symbol) # test de variantes du réseau de neurones convolutif et d’autres paramètres d’apprentissage load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_non_recouvrant.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO25.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO50.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO10_BN.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO25_BN.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO50_BN.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_BN.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO25_BN_adadelta.RData") load(file="D:/Data Mining/Logiciel R/Programmes R/Big Data/MNIST/logger_pool_recouvrant_DO25_BN_adam.RData") library(RColorBrewer) CL <- brewer.pal(10, "Spectral") library(ggplot2) trait <- 1 ggplot() + geom_line(aes(seq(2,50), logger1$eval[2:50]), colour=CL[1], size=trait, linetype=1) + geom_line(aes(seq(2,50), logger2$eval[2:50]), colour=CL[2], size=trait, linetype=2) + geom_line(aes(seq(2,50), logger3$eval[2:50]), colour=CL[3], size=trait, linetype=3) + geom_line(aes(seq(2,50), logger4$eval[2:50]), colour=CL[4], size=trait, linetype=4) + geom_line(aes(seq(2,50), logger5$eval[2:50]), colour=CL[5], size=trait, linetype=5) + geom_line(aes(seq(2,50), logger6$eval[2:50]), colour=CL[6], size=trait, linetype=6) + geom_line(aes(seq(2,50), logger7$eval[2:50]), colour=CL[7], size=trait, linetype=7) + geom_line(aes(seq(2,50), logger10$eval[2:50]), colour=CL[8], size=trait, linetype=8) + geom_line(aes(seq(2,50), logger11$eval[2:50]), colour=CL[9], size=trait, linetype=9) + geom_line(aes(seq(2,50), logger12$eval[2:50]), colour=CL[10], size=trait, linetype=10) + xlab("# passages") + ylab("taux de bon classement") + theme_bw() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 18)) + annotate("segment", x=15, xend=18, y=.971, yend=.971, colour=CL[1], size=trait, linetype=1) + annotate("text", x=20, hjust="left", y=.971, label="Pool non recouvrant", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.972, yend=.972, colour=CL[2], size=trait, linetype=2) + annotate("text", x=20, hjust="left", y=.972, label="Pool recouvrant", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.973, yend=.973, colour=CL[3], size=trait, linetype=3) + annotate("text", x=20, hjust="left", y=.973, label="Pool recouvrant DO50", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.974, yend=.974, colour=CL[4], size=trait, linetype=4) + annotate("text", x=20, hjust="left", y=.974, label="Pool recouvrant DO25", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.975, yend=.975, colour=CL[5], size=trait, linetype=6) + annotate("text", x=20, hjust="left", y=.975, label="Pool recouvrant DO50 BN", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.976, yend=.976, colour=CL[6], size=trait, linetype=6) + annotate("text", x=20, hjust="left", y=.976, label="Pool recouvrant DO25 BN", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.977, yend=.977, colour=CL[7], size=trait, linetype=7) + annotate("text", x=20, hjust="left", y=.977, label="Pool recouvrant BN", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.978, yend=.978, colour=CL[8], size=trait, linetype=8) + annotate("text", x=20, hjust="left", y=.978, label="Pool recouvrant DO10 BN", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.979, yend=.979, colour=CL[9], size=trait, linetype=9) + annotate("text", x=20, hjust="left", y=.979, label="Pool recouvrant DO25 BN ADADELTA", colour="black", size=6) + annotate("segment", x=15, xend=18, y=.980, yend=.980, colour=CL[10], size=trait, linetype=10) + annotate("text", x=20, hjust="left", y=.980, label="Pool recouvrant DO25 BN ADAM", colour="black", size=6) # Section 6.2.3 : Reconnaissance des images CIFAR-10 avec MXNet # --------------------------------------------------------------------------------------------------------- # chargement du package keras qui contient les données CIFAR-10 # voir section 5.18 pour l'installation de Keras library(keras) # chargement du jeu de données de Keras cifar10 <- dataset_cifar10() #save(cifar10, file="D:/Data/Images/cifar10.RData") x_train <- cifar10$train$x/255 x_test <- cifar10$test$x/255 y_train <- to_categorical(cifar10$train$y, num_classes = 10) y_test <- to_categorical(cifar10$test$y, num_classes = 10) # export en csv library(data.table) library(plyr) dim(x_train) cifar <- t(adply(x_train[,,,1], c(2,3), .id=NULL)) cifar <- data.frame("label"=cifar10$train$y, cifar) head(cifar) fwrite(cifar, "D:/Data/Images/CIFAR10/cifar_data_train.csv") cifar <- t(adply(x_train[,,,2], c(2,3), .id=NULL)) cifar <- data.frame("label"=cifar10$train$y, cifar) fwrite(cifar, "D:/Data/Images/CIFAR10/green_data_train.csv") cifar <- t(adply(x_train[,,,3], c(2,3), .id=NULL)) cifar <- data.frame("label"=cifar10$train$y, cifar) fwrite(cifar, "D:/Data/Images/CIFAR10/blue_data_train.csv") # variante plus rapide avec data.table DT <- as.data.table(x_train[,,,1]) head(DT) cifar <- dcast(DT, V1 ~ V3 + V2, value.var = "value")[,-1] cifar <- data.table("label"=cifar10$train$y, cifar) dim(x_test) cifar <- t(adply(x_test[,,,1], c(2,3), .id=NULL)) cifar <- data.frame("label"=cifar10$test$y, cifar) fwrite(cifar, "D:/Data/Images/CIFAR10/cifar_data_test.csv") cifar <- t(adply(x_test[,,,2], c(2,3), .id=NULL)) cifar <- data.frame("label"=cifar10$test$y, cifar) fwrite(cifar, "D:/Data/Images/CIFAR10/green_data_test.csv") cifar <- t(adply(x_test[,,,3], c(2,3), .id=NULL)) cifar <- data.frame("label"=cifar10$test$y, cifar) fwrite(cifar, "D:/Data/Images/CIFAR10/blue_data_test.csv") # variante plus rapide avec data.table DT <- as.data.table(x_test[,,,1]) head(DT) cifar <- dcast(DT, V1 ~ V3 + V2, value.var = "value")[,-1] cifar <- data.table("label"=cifar10$test$y, cifar) # chargement des données en couleurs # la 1ère variable est "label", suivie de 1024 pixels library(data.table) # échantillon d'apprentissage cifar_red <- fread("D:/Data/Images/CIFAR10/red_data_train.csv", header=TRUE) cifar_green <- fread("D:/Data/Images/CIFAR10/green_data_train.csv", header=TRUE) cifar_blue <- fread("D:/Data/Images/CIFAR10/blue_data_train.csv", header=TRUE) # transformation des matrices en arrays (nécessaire pour l'opérateur de convolution de mxnet) : table(cifar_red[,1]) train.y <- cifar_red[1:50000,1] # le label (de 0 à 9) est le même dans les trois fichiers de couleurs train.array <- array(data = NA, dim = c(32, 32, 3, 50000)) train.array[,,1,] <- t(cifar_red[1:50000,-1]) train.array[,,2,] <- t(cifar_green[1:50000,-1]) train.array[,,3,] <- t(cifar_blue[1:50000,-1]) # échantillon de test cifar_red <- fread("D:/Data/Images/CIFAR10/red_data_test.csv", header=TRUE) cifar_green <- fread("D:/Data/Images/CIFAR10/green_data_test.csv", header=TRUE) cifar_blue <- fread("D:/Data/Images/CIFAR10/blue_data_test.csv", header=TRUE) # transformation des matrices en arrays (nécessaire pour l'opérateur de convolution de mxnet) : test.y <- cifar_red[1:10000,1] # le label (de 0 à 9) est le même dans les trois fichiers de couleurs test.array <- array(data = NA, dim = c(32, 32, 3, 10000)) test.array[,,1,] <- t(cifar_red[1:10000,-1]) test.array[,,2,] <- t(cifar_green[1:10000,-1]) test.array[,,3,] <- t(cifar_blue[1:10000,-1]) # affichage des 100 premières images (en couleur) old <- par(no.readonly = TRUE) par(mfrow = c(10,10), mai = c(0,0,0,0)) for(i in 1:100) { plot(as.raster(test.array[,,,i])) # CIFAR en couleur #plot(as.raster(0.5+train.array[,,,i])) # CIFAR en gris text(3, 0, test.y[i], cex = 2, col = gray(0.75), font=3, pos = c(3,4)) } par(old) # 6 couches de convolution avec pooling non recouvrant et dropout, et normalisation batch + 2 couches complètement connectées library(mxnet) library(dplyr) cnn8 <- mx.symbol.Variable('data') %>% # première couche de convolution mx.symbol.Convolution(kernel=c(3,3), pad=c(1,1), num_filter=64) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # deuxième couche de convolution mx.symbol.Convolution(kernel=c(3,3), pad=c(1,1), num_filter=64) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(2,2), stride=c(2,2)) %>% mx.symbol.Dropout(p=0.2) %>% # troisième couche de convolution mx.symbol.Convolution(kernel=c(3,3), pad=c(1,1), num_filter=128) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # quatrième couche de convolution mx.symbol.Convolution(kernel=c(3,3), pad=c(1,1), num_filter=128) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(2,2), stride=c(2,2)) %>% mx.symbol.Dropout(p=0.2) %>% # cinquième couche de convolution mx.symbol.Convolution(kernel=c(3,3), pad=c(1,1), num_filter=256) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # sixième couche de convolution mx.symbol.Convolution(kernel=c(3,3), pad=c(1,1), num_filter=256) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% mx.symbol.Pooling(pool_type="max", kernel=c(2,2), stride=c(2,2)) %>% mx.symbol.Dropout(p=0.2) %>% # première couche complètement connectée mx.symbol.Flatten() %>% mx.symbol.Dropout(p=0.5) %>% mx.symbol.FullyConnected(num_hidden=512) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # deuxième couche complètement connectée mx.symbol.Dropout(p=0.5) %>% mx.symbol.FullyConnected(num_hidden=512) %>% mx.symbol.BatchNorm() %>% mx.symbol.Activation(act_type="relu") %>% # couche de sortie mx.symbol.FullyConnected(num_hidden=10) %>% mx.symbol.SoftmaxOutput() # entraînement du modèle #device <- mx.cpu() # entraînement du modèle sur une seule CPU device <- lapply(0:3, function(i) { mx.cpu(i) }) # avec plusieurs CPU, modifier le paramètre kvstore ci-dessous (voir https://github.com/apache/incubator-mxnet/issues/5296) logger <- mx.metric.logger$new() mx.set.seed(235) system.time(model <- mx.model.FeedForward.create(symbol=cnn8, X=train.array, y=as.numeric(as.matrix(train.y)), ctx=device, num.round=50, kvstore = "local_allreduce_device", # paramètre kvstore = "device" par défaut mais à remplacer par "local_allreduce_device" pour utiliser plusieurs CPU ou GPU # à noter qu'avec la version 3.4.2 de R, en février 2018, il n'était pas nécessaire de modifier le paramètre kvstore, qui prenait automatiquement la valeur "local_update_cpu" array.batch.size=128, optimizer="rmsprop", eval.data=list(data=test.array, label=as.numeric(as.matrix(test.y))), wd=0.001, eval.metric = mx.metric.accuracy, epoch.end.callback = mx.callback.log.train.metric(100, logger))) # prédiction preds <- predict(model, test.array) pred.label <- max.col(t(preds)) - 1 table(as.matrix(test.y), pred.label) sum(pred.label != as.matrix(test.y)) / dim(test.y)[1] # évolution du taux de bon classement plot(logger$eval, type="l", ylim=c(0.4,1), ylab="taux de bon classement", xlab="passages") abline(h=logger$eval[which.max(logger$eval)], lty=2, col="grey") lines(logger$train, lty=3) # Section 6.2.4 : Gestion des modèles avec MXNet # --------------------------------------------------------------------------------------------------------- # sauvegarde du modèle # mx.model.save(model, xx, nnnn) crée deux fichiers : xx-symbol.json et xx-nnnn.params mx.model.save(model, 10, 1) # import du modèle avec la commande mx.model.load(xx, iteration = nnnn) model.import <- mx.model.load("50", iteration = 1) preds <- predict(model.import, test.array) # application d'un modèle à une image tracée avec Paint et enregistrée en jpeg library(imager) im <- load.image("D:/Data/chiffres/chiffre4.jpg") dim(im) plot(im) im <- 1 - im # inversion pour avoir du blanc sur fond noir resized <- resize(im, 28, 28) resized <- grayscale(resized) # suppression des couleurs pour ne conserver qu'un seul canal resized <- resized/max(resized) # normalisation des pixels entre 0 et 1 dim(resized) plot(resized) test.array <- as.array(resized) preds <- predict(model.import, test.array) (pred.label <- max.col(t(preds)) - 1) # chiffre prédit # code sous forme de fonction pour tester un ensemble de chiffres library(imager) reconnaissance <- function(x){ im <- load.image(paste("D:/Data/chiffres/",x, sep="")) im <- 1 - im # inversion pour avoir du blanc (gray(1)) sur fond noir (gray(0)) resized <- resize(im, 28, 28) resized <- grayscale(resized) # suppression des couleurs pour ne conserver qu'un seul canal resized <- resized/max(resized) # normalisation des pixels entre 0 et 1 test.array <- as.array(resized) preds <- predict(model.import, test.array) (pred.label <- max.col(t(preds)) - 1) # chiffre prédit } reconnaissance("chiffre4.jpg") sapply(0:9, function(x) {reconnaissance(paste("chiffre",x,".jpg", sep=""))}) # Section 6.3.2 : Application de Keras aux images MNIST # --------------------------------------------------------------------------------------------------------- # chargement du package keras # voir section 5.18 pour l'installation de Keras library(keras) # chargement du jeu de données c(c(x_train, y_train), c(x_test, y_test)) %<-% dataset_mnist() # dimensions des images en entrée img_rows <- 28 img_cols <- 28 num_classes <- 10 # conversion d'un tenseur 3D (ici 60000x28x28 pour l'apprentissage) en une matrice 2D pour un perceptron (MLP) x_train <- array(x_train, dim = c(nrow(x_train), img_rows*img_cols)) x_test <- array(x_test, dim = c(nrow(x_test), img_rows*img_cols)) # conversion d'un tenseur 3D (ici 60000x28x28 pour l'apprentissage) en un tenseur 4D pour un CNN x_train <- array(x_train, c(nrow(x_train), img_rows, img_cols, 1)) x_test <- array(x_test, c(nrow(x_test), img_rows, img_cols, 1)) input_shape <- c(img_rows, img_cols, 1) # les valeurs RGB sont normalisées dans [0,1] x_train <- x_train / 255 x_test <- x_test / 255 # conversion des numéros de classes en matrices binaires (une colonne par classe) y_train <- to_categorical(y_train, num_classes) y_test <- to_categorical(y_test, num_classes) # modèle séquentiel keras model <- keras_model_sequential() # MLP avec une couche en entrée [784 unités], 1 couche cachée [100 unités] avec un dropout de 0.5, et une couche de sortie [10 unités] model %>% layer_dense(units = 100, input_shape = 784) %>% layer_dropout(rate=0.5) %>% layer_activation(activation = 'relu') %>% layer_dense(units = 10) %>% layer_activation(activation = 'softmax') # autre modèle (MLP à deux couches cachées) model %>% layer_dense(units = 256, activation = 'relu', input_shape = c(784)) %>% layer_dropout(rate = 0.4) %>% layer_dense(units = 128, activation = 'relu') %>% layer_dropout(rate = 0.3) %>% layer_dense(units = 10, activation = 'softmax') # autre modèle (CNN) model %>% layer_conv_2d(filter = 10, kernel_size = c(5,5), input_shape = input_shape) %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(3,3), strides = c(2,2)) %>% layer_conv_2d(filter = 20, kernel_size = c(5,5)) %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(3,3), strides = c(2,2)) %>% layer_flatten() %>% layer_dropout(0.5, seed=235) %>% layer_dense(300) %>% layer_activation("relu") %>% layer_dense(num_classes) %>% layer_activation("softmax") # autre modèle (CNN) : celui du livre model %>% layer_conv_2d(filter = 20, kernel_size = c(5,5), input_shape = input_shape) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(3,3), strides = c(2,2)) %>% layer_conv_2d(filter = 60, kernel_size = c(5,5)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(3,3), strides = c(2,2)) %>% layer_flatten() %>% layer_dropout(0.5, seed=235) %>% layer_dense(600) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_dense(num_classes) %>% layer_activation("softmax") summary(model) # compilation du modèle model %>% compile(loss = 'categorical_crossentropy', optimizer = 'sgd', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = 'adam', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = 'rmsprop', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = 'adadelta', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = 'adagrad', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = optimizer_rmsprop(lr = 0.0001, decay = 1e-6), metrics = 'accuracy') model %>% compile(loss = 'categorical_crossentropy', optimizer = optimizer_adadelta(decay = 1e-5), metrics = 'accuracy') # autre fonction de perte (pour une régression) : mse # apprentissage du modèle sur l'échantillon d'apprentissage system.time(history <- fit(model, x_train, y_train, epochs = 50, batch_size = 128, shuffle = TRUE, validation_data = list(x_test, y_test), verbose=0)) plot(history) # évaluation du modèle sur l'échantillon de validation model %>% evaluate(x_test, y_test, verbose = 0) # prédictions #pred <- predict_proba(model, x_test) pred <- predict_classes(model, x_test) test.label <- max.col(y_test) - 1 table(test.label, pred) mean(test.label != pred) # Section 6.3.3 : Application de modèles préentraînés # --------------------------------------------------------------------------------------------------------- # enregistrement du modèle save_model_hdf5(model, "full_model.h5") # modèle entier enregistré comme fichier binaire save_model_weights_hdf5(model, "weights_model.h5") # ensemble des poids du modèle enregistré comme fichier binaire # lecture du modèle model <- load_model_hdf5("full_model.h5") load_model_weights_hdf5(model, "weights_model.h5") # application d'un modèle à une image tracée avec Paint et enregistrée en jpeg library(imager) im <- load.image("D:/Data/chiffres/chiffre4.jpg") dim(im) im <- 1 - im # inversion pour avoir du blanc (gray(1)) sur fond noir (gray(0)) resized <- resize(im, 28, 28) resized <- grayscale(resized) # suppression des couleurs pour ne conserver qu'un seul canal resized <- resized/max(resized) # normalisation des pixels entre 0 et 1 dim(resized) test.array <- array(resized, c(1,28,28,1)) test.array test.array[1,,,] # chiffre prédit model <- load_model_hdf5("full_model.h5") # permutation hauteur et la largeur de l'image dans le tenseur (pred <- predict_classes(model, aperm(test.array,c(1,3,2,4)))) # code sous forme de fonction pour tester un ensemble de chiffres (avec fonctions de preprocessing de imager) reconnaissance <- function(x){ im <- load.image(paste("D:/Data/chiffres/",x, sep="")) im <- 1 - im # inversion pour avoir du blanc (gray(1)) sur fond noir (gray(0)) resized <- resize(im, 28, 28) resized <- grayscale(resized) # suppression des couleurs pour ne conserver qu'un seul canal resized <- resized/max(resized) # normalisation des pixels entre 0 et 1 test.array <- array(resized, c(1,28,28,1)) # permutation hauteur et la largeur de l'image dans le tenseur (pred <- predict_classes(model, aperm(test.array,c(1,3,2,4)))) #(pred <- predict_classes(model, test.array)) } #reconnaissance("chiffre1.jpg") sapply(0:9, function(x) {reconnaissance(paste("chiffre",x,".jpg", sep=""))}) # code sous forme de fonction pour tester un ensemble de chiffres (avec fonctions de preprocessing de keras) reconnaissance <- function(x){ img <- image_load(paste("D:/Data/chiffres/",x, sep=""), target_size = c(28,28), grayscale = TRUE) x <- image_to_array(img) # conversion de l'image en tenseur 3D x <- array_reshape(1-(x/255), c(1, dim(x))) # normalisation des pixels entre 0 et 1 et inversion blanc et noir # tenseur 4D avec ajout 1ère composante = 1 car une seule image preds <- model %>% predict(x) which(preds==max(preds))-1 } #reconnaissance("chiffre4.jpg") sapply(0:9, function(x) {reconnaissance(paste("chiffre",x,".jpg", sep=""))}) # chargement d'un modèle préentraîné model_vgg <- application_vgg16(weights = "imagenet", include_top = TRUE) model_vgg <- application_vgg19(weights = "imagenet", include_top = TRUE) model_res <- application_resnet50(weights = "imagenet", include_top = TRUE) model_res <- application_inception_v3(weights = "imagenet", include_top = TRUE) model_res <- application_inception_resnet_v2(weights = "imagenet", include_top = TRUE) model_res <- application_xception(weights = "imagenet", include_top = TRUE) summary(model_vgg) img_path <- "D:/Data/Images/berger_allemand.jpg" img <- image_load(img_path, target_size = c(224,224)) # chargement de l'image x <- image_to_array(img) # conversion de l'image en tenseur 3D x <- array_reshape(x, c(1, dim(x))) # tenseur 4D avec ajout 1ère composante = 1 car une seule image # même préprocessing que pour le concours ImageNet : soustraction de la moyenne RGB pour chaque pixel x <- imagenet_preprocess_input(x) preds <- model_vgg %>% predict(x) # calcul de probabilités pour chacune des 1000 classes dim(preds) imagenet_decode_predictions(preds, top = 3)[[1]] # transfer learning : extraction de caractéristiques base_model <- application_vgg16(weights = "imagenet", include_top = FALSE, input_shape = c(32, 32, 3)) features <- base_model %>% predict(x) dim(features) base_model # ajout d'une couche complètement connectée sur la base convolutive model <- keras_model_sequential() %>% base_model %>% layer_flatten() %>% #layer_global_average_pooling_2d() %>% layer_dense(units = 128, activation = "relu") %>% layer_dropout(rate = 0.5) %>% #layer_dense(units = 1, activation = "sigmoid") # prédiction binaire (type chat / chien) layer_dense(units = 10, activation = "softmax") # prédiction multiclasses (type MNIST ou CIFAR) model freeze_weights(base_model) unfreeze_weights(base_model, from = "block3_conv1") model # compilation pour prédiction binaire #model %>% compile(loss = 'binary_crossentropy', optimizer = optimizer_rmsprop(lr = 0.001), metrics = 'accuracy') # compilation pour prédiction multiclasses model %>% compile(loss = 'categorical_crossentropy', optimizer = optimizer_rmsprop(lr = 0.001, decay = 1e-3), metrics = 'accuracy') # apprentissage sur échantillon system.time(history <- fit(model, x_train, y_train, epochs = 30, batch_size = 128, shuffle = TRUE, validation_data = list(x_test, y_test), verbose=1)) # apprentissage sur lots d'images générés #system.time(history <- model %>% fit_generator(train_generator, steps_per_epoch = 80, epochs = 30, validation_data = test_generator, validation_steps = 40)) plot(history) # Section 6.3.4 : Application de Keras aux images CIFAR-10 # --------------------------------------------------------------------------------------------------------- # chargement du package keras # voir section 5.18 pour l'installation de Keras library(keras) # chargement du jeu de données cifar10 <- dataset_cifar10() # facultatif : normalisation for(i in 1:3){ mea <- mean(cifar10$train$x[,,,i]) sds <- sd(cifar10$train$x[,,,i]) cifar10$train$x[,,,i] <- (cifar10$train$x[,,,i] - mea) / sds cifar10$test$x[,,,i] <- (cifar10$test$x[,,,i] - mea) / sds } x_train <- cifar10$train$x x_test <- cifar10$test$x # avec normalisation basique x_train <- cifar10$train$x/255 x_test <- cifar10$test$x/255 y_train <- to_categorical(cifar10$train$y, num_classes = 10) y_test <- to_categorical(cifar10$test$y, num_classes = 10) # affichage des 100 premières images old <- par(no.readonly = TRUE) par(mfrow = c(10,10), mai = c(0,0,0,0)) for(i in 1:100) { plot(as.raster(x_train[i,,,])) text(3, 0, cifar10$train$y[i], cex = 2, col = gray(0.75), font=3, pos = c(3,4)) } par(old) # modèle séquentiel keras à 4 couches de convolution (code absent du livre) model <- keras_model_sequential() %>% # première couche de convolution, alimentée par des images de 32x32 pixels layer_conv_2d(filter = 32, kernel_size = c(3,3), padding = "same", input_shape = c(32, 32, 3)) %>% # layer_batch_normalization() %>% layer_activation("relu") %>% # deuxième couche de convolution layer_conv_2d(filter = 32, kernel_size = c(3,3)) %>% # layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(2,2)) %>% layer_dropout(0.2) %>% # 3e couche de convolution layer_conv_2d(filter = 32, kernel_size = c(3,3), padding = "same") %>% # layer_batch_normalization() %>% layer_activation("relu") %>% # 4e couche de convolution layer_conv_2d(filter = 32, kernel_size = c(3,3)) %>% # layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(2,2)) %>% layer_dropout(0.2) %>% # première couche complètement connectée layer_flatten() %>% layer_dropout(0.5) %>% layer_dense(512) %>% # layer_batch_normalization() %>% layer_activation("relu") %>% # couche de sortie layer_dense(10) %>% layer_activation("softmax") summary(model) # modèle séquentiel keras à 6 couches de convolution (p. 198) model <- keras_model_sequential() %>% # première couche de convolution layer_conv_2d(filter = 64, kernel_size = c(3,3), padding = "same", input_shape = c(32, 32, 3)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% # deuxième couche de convolution layer_conv_2d(filter = 64, kernel_size = c(3,3)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(2,2)) %>% layer_dropout(0.2) %>% # 3e couche de convolution layer_conv_2d(filter = 128, kernel_size = c(3,3), padding = "same") %>% layer_batch_normalization() %>% layer_activation("relu") %>% # 4e couche de convolution layer_conv_2d(filter = 128, kernel_size = c(3,3)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(2,2)) %>% layer_dropout(0.2) %>% # 5e couche de convolution layer_conv_2d(filter = 256, kernel_size = c(3,3), padding = "same") %>% layer_batch_normalization() %>% layer_activation("relu") %>% # 6e couche de convolution layer_conv_2d(filter = 256, kernel_size = c(3,3)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_2d(pool_size = c(2,2)) %>% layer_dropout(0.2) %>% # première couche complètement connectée layer_flatten() %>% layer_dropout(0.5) %>% layer_dense(units = 512) %>% layer_batch_normalization() %>% layer_activation("relu") %>% # deuxième couche complètement connectée layer_dropout(0.5) %>% layer_dense(units = 512) %>% layer_batch_normalization() %>% layer_activation("relu") %>% # couche de sortie layer_dense(10) %>% # sortie layer_activation("softmax") summary(model) # compilation du modèle model %>% compile(loss = 'categorical_crossentropy', optimizer = optimizer_rmsprop(lr = 1e-3, decay = 1e-3), metrics = 'accuracy') # apprentissage du modèle sur l'échantillon d'apprentissage (sans data augmentation) system.time(history <- fit(model, x_train, y_train, epochs = 50, batch_size = 128, shuffle = TRUE, validation_data = list(x_test, y_test), verbose=1)) plot(history) # évaluation du modèle sur l'échantillon de validation model %>% evaluate(x_test, y_test, verbose = 0) # prédictions pred <- predict_classes(model, x_test) # classes prédites #pred <- predict_proba(model, x_test) # proba de chacune des classes prédites (une par colonne) test.label <- max.col(y_test) - 1 table(test.label, pred) mean(test.label != pred) # affichage des 100 premières erreurs de classement err100 <- head(which(test.label != pred), 100) old <- par(no.readonly = TRUE) par(mfrow = c(10,10), mai = c(0,0,0,0)) for(i in err100) { plot(as.raster(x_test[i,,,])) text(3, 0, test.label[i], cex = 2, col = gray(0.75), font=3, pos = c(3,4)) text(3, 20, pred[i], cex = 2, col = gray(0.25), pos = c(3,4)) } par(old) # augmentation de données (p. 204) # paramètres pour générer des observations différentes sans trop les déformer datagen <- image_data_generator( rotation_range = 15, width_shift_range = 0.15, height_shift_range = 0.15, shear_range = 0.15, zoom_range = 0.15, horizontal_flip = FALSE ) # apprentissage du modèle sur l'échantillon d'apprentissage (avec data augmentation) system.time(history <- model %>% fit_generator( flow_images_from_data(x_train, y_train, generator = datagen, batch_size = 128), steps_per_epoch = as.integer(nrow(x_train)/128), epochs = 50, validation_data = list(x_test, y_test) # , callbacks = callback_early_stopping(monitor = 'val_loss', patience = 3) )) # évaluation du modèle sur l'échantillon de validation model %>% evaluate(x_test, y_test, verbose = 0) # chargement d'un modèle préentraîné (code absent du livre) base_model <- application_vgg16(include_top = FALSE, weights = "imagenet", input_shape = c(32, 32, 3)) base_model <- application_mobilenet(include_top = FALSE, weights = "imagenet", input_shape = c(128, 128, 3)) base_model <- application_resnet50(include_top = FALSE, weights = "imagenet", input_shape = c(224, 224, 3)) base_model # complétion du modèle avec deux couches complètement connectées choisies identiques à celles des modèles précédents model <- keras_model_sequential() %>% base_model %>% layer_flatten() %>% # première couche complètement connectée layer_dropout(0.5) %>% layer_dense(units = 512) %>% layer_batch_normalization() %>% layer_activation("relu") %>% # deuxième couche complètement connectée layer_dropout(0.5) %>% layer_dense(units = 512) %>% layer_batch_normalization() %>% layer_activation("relu") %>% # couche de sortie layer_dense(units = 10, activation = "softmax") model # complétion du modèle avec une couche complètement connectée model <- keras_model_sequential() %>% base_model %>% layer_flatten() %>% layer_dropout(0.5) %>% layer_dense(units = 512, kernel_regularizer = regularizer_l2(0.0001)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_dense(units = 10, activation = "softmax") model # autre complétion du modèle avec une couche complètement connectée # cas d'images à traiter plus petites que l'input_shape du modèle pré-entraîné model <- keras_model_sequential() %>% layer_upsampling_2d %>% # pour passer d'images m x m à des images 2m x 2m en répétant les images et les colonnes layer_upsampling_2d %>% layer_upsampling_2d %>% base_model %>% layer_flatten() %>% layer_batch_normalization() %>% layer_dense(units = 128, activation = "relu") %>% layer_dropout(0.5) %>% #layer_batch_normalization() %>% #layer_dense(units = 64, activation = "relu") %>% #layer_dropout(0.5) %>% layer_batch_normalization() %>% layer_dense(units = 10, activation = "softmax") model # modèle en partant d'une couche intermédiaire de VGG16 predictions <- get_layer(base_model, 'block3_pool')$output %>% layer_flatten() %>% layer_dropout(0.5) %>% layer_dense(units = 512, kernel_regularizer = regularizer_l2(0.0001)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_dense(units = 10, activation = "softmax") model <- keras_model(inputs = base_model$input, outputs = predictions) model # équivalent à ce qui précède : modèle en partant d'une couche intermédiaire de VGG16 intermediate_layer_model <- keras_model(inputs = base_model$input, outputs = get_layer(base_model, 'block3_pool')$output) model <- keras_model_sequential() %>% intermediate_layer_model %>% layer_flatten() %>% layer_dropout(0.5) %>% layer_dense(units = 512, kernel_regularizer = regularizer_l2(0.0001)) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_dense(units = 10, activation = "softmax") model # gel des couches convolutives freeze_weights(base_model) model # équivalent for (layer in base_model$layers){ layer$trainable = FALSE } # après un 1er apprentissage des couches denses, on peut dégeler une partie des couches convolutives unfreeze_weights(base_model, from = "block5_conv1") unfreeze_weights(base_model, from = "block4_conv1") unfreeze_weights(base_model, from = "block3_pool") model # compilation du modèle model %>% compile(loss = 'categorical_crossentropy', optimizer = optimizer_rmsprop(lr = 1e-3, decay = 1e-3), metrics = 'accuracy') model %>% compile(loss = 'binary_crossentropy', optimizer = optimizer_rmsprop(lr = 2e-5), metrics = 'accuracy') # apprentissage du modèle sans data augmentation system.time(history <- fit(model, x_train, y_train, epochs = 10, batch_size = 128, shuffle = TRUE, validation_data = list(x_test, y_test), verbose=1)) # apprentissage du modèle avec data augmentation system.time(history <- model %>% fit_generator( flow_images_from_data(x_train, y_train, generator = datagen, shuffle = TRUE, batch_size = 128), steps_per_epoch = as.integer(nrow(x_train)/128), epochs = 50, validation_data = list(x_test, y_test) )) plot(history) model %>% evaluate(x_test, y_test, verbose = 0) # prédictions pred <- predict_classes(model, x_test) test.label <- max.col(y_test) - 1 table(test.label, pred) mean(test.label != pred) # Section 6.3.5 : Application de Keras au classement de photos de chats et de chiens # --------------------------------------------------------------------------------------------------------- # Lecture des données # -------------------------------------------------------------------- # # On peut sauter l'étape suivante si les images sont déjà rangées dans les bons répertoires : # un pour les images d’apprentissage et un pour les images de test, chacun contenant deux sous-répertoires : pour les chiens et pour les chats original_dataset_dir <- "D:/Data/DogsCats" base_dir <- "D:/Data/DogsCats/miniDogsCats" dir.create(base_dir) train_dir <- file.path(base_dir, "train") dir.create(train_dir) test_dir <- file.path(base_dir, "test") dir.create(test_dir) train_cats_dir <- file.path(train_dir, "cats") dir.create(train_cats_dir) train_dogs_dir <- file.path(train_dir, "dogs") dir.create(train_dogs_dir) test_cats_dir <- file.path(test_dir, "cats") dir.create(test_cats_dir) test_dogs_dir <- file.path(test_dir, "dogs") dir.create(test_dogs_dir) fnames <- paste0("cat.", 1:2500, ".jpg") file.copy(file.path(original_dataset_dir, "train", fnames), file.path(train_cats_dir)) fnames <- paste0("cat.", 3001:4000, ".jpg") file.copy(file.path(original_dataset_dir, "train", fnames), file.path(test_cats_dir)) fnames <- paste0("dog.", 1:2500, ".jpg") file.copy(file.path(original_dataset_dir, "train", fnames), file.path(train_dogs_dir)) fnames <- paste0("dog.", 3001:4000, ".jpg") file.copy(file.path(original_dataset_dir, "train", fnames), file.path(test_dogs_dir)) # Modèle simple (p. 207) # -------------------------------------------------------------------- # img_rows = 128 img_cols = 128 input_shape <- c(img_rows, img_cols, 3) model <- keras_model_sequential() %>% layer_conv_2d(filters = 32, kernel_size = c(3,3), activation = 'relu', input_shape = input_shape) %>% 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 = 128, kernel_size = c(3,3), activation = 'relu') %>% layer_max_pooling_2d(pool_size = c(2, 2)) %>% layer_flatten() %>% layer_dense(units = 128, activation = 'relu') %>% layer_dropout(rate = 0.5) %>% layer_dense(units = 1, activation = 'sigmoid') summary(model) # compilation du modèle model %>% compile( loss = 'binary_crossentropy', #loss = 'categorical_crossentropy', #optimizer = optimizer_adadelta(), optimizer = optimizer_rmsprop(lr=0.001), #optimizer = optimizer_adam(lr=0.001, beta_1=0.9, beta_2=0.9, decay=0.1), metrics = c('accuracy') ) train_datagen <- image_data_generator(rescale = 1/255) test_datagen <- image_data_generator(rescale = 1/255) base_dir <- "D:/Data/DogsCats/miniDogsCats" train_dir <- file.path(base_dir, "train") test_dir <- file.path(base_dir, "test") train_generator <- flow_images_from_directory(train_dir, train_datagen, target_size = c(img_rows,img_cols), shuffle = TRUE, seed = 235, batch_size = 64, class_mode = "binary") # option shuffle = FALSE sur les données de test en vue du calcul plus bas de la matrice de confusion test_generator <- flow_images_from_directory(test_dir, test_datagen, target_size = c(img_rows,img_cols), shuffle = FALSE, batch_size = 64, class_mode = "binary") # entraînement du modèle system.time(history <- model %>% fit_generator(train_generator, steps_per_epoch = 80, epochs = 30, validation_data = test_generator, validation_steps = 30)) plot(history) # évaluation du modèle sur l'échantillon de validation model %>% evaluate_generator(test_generator, steps = 30) # prédictions pred <- predict_generator(model, test_generator, steps = 30) # probas prédites table(head(test_generator$classes,length(pred)), (pred>0.5)) # Modèle avec data augmentation (p. 209) # -------------------------------------------------------------------- # datagen <- image_data_generator( rescale = 1/255, rotation_range = 15, width_shift_range = 0.15, height_shift_range = 0.15, shear_range = 0.15, zoom_range = 0.15, horizontal_flip = FALSE ) # apprentissage du modèle sur l'échantillon d'apprentissage (avec data augmentation) train_generator <- flow_images_from_directory(train_dir, datagen, target_size = c(img_rows,img_cols), shuffle = TRUE, seed = 235, batch_size = 64, class_mode = "binary") system.time(history <- model %>% fit_generator(train_generator, steps_per_epoch = 80, epochs = 30, validation_data = test_generator, validation_steps = 30)) plot(history) # Modèle avec modèle préentraîné (p. 211) # -------------------------------------------------------------------- # # chargement d'un modèle préentraîné base_model <- application_vgg16(include_top = FALSE, weights = "imagenet", input_shape = input_shape) base_model <- application_inception_v3(include_top = FALSE, weights = 'imagenet', input_shape = input_shape) base_model <- application_inception_resnet_v2(include_top = FALSE, weights = 'imagenet', input_shape = input_shape) base_model <- application_xception(include_top = FALSE, weights = 'imagenet', input_shape = input_shape) base_model <- application_mobilenet(include_top = FALSE, weights = 'imagenet', input_shape = input_shape) base_model # ajout d'une couche complètement connectée sur la base convolutive model <- keras_model_sequential() %>% base_model %>% layer_flatten() %>% #layer_global_average_pooling_2d() %>% layer_dense(units = 128, activation = "relu") %>% layer_dropout(rate = 0.5) %>% layer_dense(units = 1, activation = "sigmoid") model freeze_weights(base_model) #unfreeze_weights(base_model, from = "block5_conv1") model model %>% compile(loss = 'binary_crossentropy', optimizer = optimizer_rmsprop(lr = 0.001), metrics = 'accuracy') system.time(history <- model %>% fit_generator(train_generator, steps_per_epoch = 80, epochs = 30, validation_data = test_generator, validation_steps = 30)) plot(history) # Section 6.6 : PyTorch # --------------------------------------------------------------------------------------------------------- # le code qui suit est en Python # chargement des packages import torch from torchvision import datasets, transforms # paramètres généraux batch_size = 64 #batch_size_test = 1000 # chargement des données FashionMNIST # transformation des données en entrée en tenseur PyTorch transfo = transforms.Compose([transforms.ToTensor()]) # paramètres de normalisation des données en entrée #transfo = transforms.Compose([transforms.ToTensor(), transforms.Normalize((0.1307,), (0.3081,))]) # téléchargement et chargement des données d'apprentissage # on indique la source des données, les transformations à leur faire subir, et le fait qu'il faut les télécharger # ensuite le data loader fournit les données au réseau, par lots de la taille indiquée, éventuellement en les mélangeant train_dataset = datasets.FashionMNIST('/files/', train=True, transform=transfo, download=True) train_loader = torch.utils.data.DataLoader(dataset=train_dataset, batch_size=batch_size, shuffle=True) # téléchargement et chargement des données de test test_dataset = datasets.FashionMNIST('/files/', train=False, transform=transfo, download=True) test_loader = torch.utils.data.DataLoader(dataset=test_dataset, batch_size=batch_size, shuffle=False) # affichage d'images examples = enumerate(train_loader) batch_idx, (example_data, example_targets) = next(examples) example_data.shape import matplotlib.pyplot as plt labels_dict = {0 : 'T-Shirt', 1 : 'Trouser', 2 : 'Pullover', 3 : 'Dress', 4 : 'Coat', 5 : 'Sandal', 6 : 'Shirt', 7 : 'Sneaker', 8 : 'Bag', 9 : 'Ankle Boot' } print(labels_dict[3]) labels_map = {0 : 'T-Shirt', 1 : 'Trouser', 2 : 'Pullover', 3 : 'Dress', 4 : 'Coat', 5 : 'Sandal', 6 : 'Shirt', 7 : 'Sneaker', 8 : 'Bag', 9 : 'Ankle Boot'}; fig = plt.figure() for i in range(6): plt.subplot(2,3,i+1) plt.tight_layout() plt.imshow(example_data[i][0], cmap='gray', interpolation='none') plt.title(labels_dict[int(example_targets[i])]) plt.xticks([]) plt.yticks([]) fig # construction du réseau import torch.nn as nn import torch.nn.functional as F import torch.optim as optim # modèle à 2 couches de convolution # le modèle est créé comme sous-classe de la classe nn.Module # la classe nn.Module est une classe de PyTorch qui contient tout ce qui est nécessaire pour créer un réseau de neurones profond # pour créer une classe, il faut définir une fonction __init__ # - qui a des caractéristiques conv1, conv2... créées à partir de la méthode Conv2d de la classe nn.Module # - et une méthode forward qui se substitue à la méthode forward de la classe nn.Module class ConvNet(nn.Module): def __init__(self): super(ConvNet, self).__init__() self.conv1 = nn.Conv2d(1, 32, kernel_size=5, stride=1, padding=0) self.conv1_bn = nn.BatchNorm2d(32) self.conv2 = nn.Conv2d(32, 64, kernel_size=5, stride=1, padding=0) self.conv2_bn = nn.BatchNorm2d(64) self.fc1 = nn.Linear(4*4*64, 256) # en sortie des couches de convolution, on a 256 cartes de dim 4 x 4 self.fc1_bn = nn.BatchNorm1d(256) # 4 x 4 vient de 28 -> 24 -> 12 -> 8 -> 4 self.fc2 = nn.Linear(256, 10) def forward(self, x): x = F.relu(self.conv1_bn(self.conv1(x))) x = F.dropout(x, p=0.1, training=self.training) x = F.max_pool2d(x, kernel_size=2, stride=2) x = F.relu(self.conv2_bn(self.conv2(x))) x = F.max_pool2d(x, kernel_size=2, stride=2) x = F.dropout(x, p=0.1, training=self.training) x = x.view(-1, 4*4*64) #x = F.relu(self.fc1(x)) x = F.relu(self.fc1_bn(self.fc1(x))) x = F.dropout(x, p=0.2, training=self.training) x = self.fc2(x) return F.log_softmax(x, dim=1) # autre modèle à 2 couches de convolution (absent du livre) class ConvNet(nn.Module): def __init__(self): super(ConvNet, self).__init__() self.conv1 = nn.Conv2d(1, 32, kernel_size=5) self.conv2 = nn.Conv2d(32, 64, kernel_size=5) self.conv2_bn = nn.BatchNorm2d(64) self.fc1 = nn.Linear(4*4*64, 256) self.fc1_bn = nn.BatchNorm1d(256) self.fc2 = nn.Linear(256, 10) def forward(self, x): x = F.relu(self.conv1(x)) x = F.dropout(x, p=0.1, training=self.training) x = F.max_pool2d(x,2) #x = F.relu(self.conv2(x)) x = F.relu(self.conv2_bn(self.conv2(x))) x = F.max_pool2d(x,2) x = F.dropout(x, p=0.1, training=self.training) x = x.view(-1,4*4*64) #x = F.relu(self.fc1(x)) x = F.relu(self.fc1_bn(self.fc1(x))) x = F.dropout(x, p=0.1, training=self.training) x = self.fc2(x) return F.log_softmax(x, dim=1) # dim=1 pour appliquer softmax par ligne # modèle à 3 couches de convolution (absent du livre) conv_units = [32,32,64] dense_units = [256] output_units = 10 class ConvNet(nn.Module): def __init__(self): super(ConvNet, self).__init__() self.conv1 = nn.Conv2d(1, conv_units[0], kernel_size=5, stride=1, padding=0) self.conv2 = nn.Conv2d(conv_units[0], conv_units[1], kernel_size=5, stride=1, padding=0) self.conv2_bn = nn.BatchNorm2d(conv_units[1]) self.conv3 = nn.Conv2d(conv_units[1], conv_units[2], kernel_size=5, stride=1, padding=0) self.conv3_bn = nn.BatchNorm2d(conv_units[2]) self.fc1 = nn.Linear(3*3*conv_units[2], dense_units[0]) # en sortie des couches de convolution, on a 256 cartes de dim 4 x 4 self.fc1_bn = nn.BatchNorm1d(dense_units[0]) # 3 x 3 vient de 28 -> 24 -> 20 -> 10 -> 6 -> 3 self.fc2 = nn.Linear(dense_units[0], output_units) def forward(self, x): x = F.relu(self.conv1(x)) #x = F.dropout(x, p=0.1, training=self.training) #x = F.relu(self.conv2(x)) x = F.relu(self.conv2_bn(self.conv2(x))) x = F.max_pool2d(x, kernel_size=2, stride=2) #x = F.dropout(x, p=0.1, training=self.training) x = F.relu(self.conv3_bn(self.conv3(x))) x = F.max_pool2d(x, kernel_size=2, stride=2) #x = F.dropout(x, p=0.1, training=self.training) x = x.view(-1, 3*3*conv_units[2]) #x = F.relu(self.fc1(x)) x = F.relu(self.fc1_bn(self.fc1(x))) x = F.dropout(x, p=0.2, training=self.training) x = self.fc2(x) return F.log_softmax(x, dim=1) # modèle à 3 couches de convolution class ConvNet(nn.Module): def __init__(self): super(ConvNet, self).__init__() self.conv1 = nn.Conv2d(1, 32, kernel_size=5, stride=1, padding=0) self.conv1_bn = nn.BatchNorm2d(32) self.conv2 = nn.Conv2d(32, 32, kernel_size=5, stride=1, padding=0) self.conv2_bn = nn.BatchNorm2d(32) self.conv3 = nn.Conv2d(32, 64, kernel_size=5, stride=1, padding=0) self.conv3_bn = nn.BatchNorm2d(64) self.fc1 = nn.Linear(3*3*64, 256) # en sortie des couches de convolution, on a 256 cartes de dim 3 x 3 self.fc1_bn = nn.BatchNorm1d(256) # 3 x 3 vient de 28 -> 24 -> 20 -> 10 -> 6 -> 3 self.fc2 = nn.Linear(256, 10) def forward(self, x): #x = F.relu(self.conv1(x)) x = F.relu(self.conv1_bn(self.conv1(x))) x = F.dropout(x, p=0.1, training=self.training) #x = F.relu(self.conv2(x)) x = F.relu(self.conv2_bn(self.conv2(x))) x = F.max_pool2d(x, kernel_size=2, stride=2) x = F.dropout(x, p=0.1, training=self.training) x = F.relu(self.conv3_bn(self.conv3(x))) x = F.max_pool2d(x, kernel_size=2, stride=2) x = F.dropout(x, p=0.1, training=self.training) x = x.flatten(1) #x = x.view(x.size(0), -1) #x = x.view(-1, 3*3*64) #x = F.relu(self.fc1_bn(x)) x = F.relu(self.fc1_bn(self.fc1(x))) x = F.dropout(x, p=0.2, training=self.training) x = self.fc2(x) return F.log_softmax(x, dim=1) # création d'une instance de la classe ConvNet model = ConvNet() #optimizer = optim.SGD(model.parameters(), lr=0.1, momentum=0.9) optimizer = torch.optim.Adam(model.parameters()) # lr=0.001, betas=(0.9,0.999)) optimizer = torch.optim.Adamax(model.parameters()) # lr=0.001, betas=(0.9,0.999)) optimizer = torch.optim.RMSprop(model.parameters()) # lr=0.001, betas=(0.9,0.999)) optimizer = torch.optim.RMSprop(model.parameters(), lr=0.001) optimizer = torch.optim.ASGD(model.parameters()) # lr=0.001, betas=(0.9,0.999)) optimizer = optim.Rprop(model.parameters()) # Entraînement et test du modèle n_epochs = 10 import time debut = time.time() loss_list = [] for epoch in range(n_epochs): model.train() correct = 0 total = 0 for images, labels in train_loader: # propagation output = model(images) loss = F.cross_entropy(output, labels) loss_list.append(loss.item()) # rétropropagation optimizer.zero_grad() # mise à 0 du gradient loss.backward() optimizer.step() # suivi de la précision _, predicted = torch.max(output.data, 1) total += labels.size(0) correct += (predicted == labels).sum().item() print('Epoch [{}/{}], Train Loss: {:.4f}, Train Accuracy: {:.2f}%, ' .format(epoch + 1, n_epochs, loss.item(), (correct / total) * 100), end='') model.eval() with torch.no_grad(): correct = 0 total = 0 for images, labels in test_loader: output = model(images) _, predicted = torch.max(output.data, 1) total += labels.size(0) correct += (predicted == labels).sum().item() print('Test Accuracy: {:.2f}%'.format((correct / total) * 100)) print('Temps de calcul : {:.0f} s.'.format(time.time() - debut)) # enregistrement du modèle MODEL_STORE_PATH = 'D:/Data Mining/Logiciel R/Programmes R/Big Data/' torch.save(model.state_dict(), MODEL_STORE_PATH + 'conv_net_model.ckpt') # fin du code Python # --------------------------------------------------------------------------------------------------------- # Chapitre 7 : La reconnaissance de l’écriture manuscrite # --------------------------------------------------------------------------------------------------------- # Section 7.3 : Récupération et mise en forme des données # --------------------------------------------------------------------------------------------------------- # les images sont de taille 28 x 28 pixels # le fichier train.csv a 785 colonnes : # - label = chiffre # - pixel0 à pixel783 sont les 784 = 28 x 28 pixels de l'image du chiffre # test.csv a 785 colonnes et 28000 lignes (chiffres) # les valeurs des pixels sont des niveaux de gris entre 0 et 255 # données directement disponibles au format csv : http://www.pjreddie.com/projects/mnist-in-csv/ # lecture des données train <- read.csv("D:/Data/MNIST/mnist_train.csv", header=TRUE) test <- read.csv("D:/Data/MNIST/mnist_test.csv", header=TRUE) # affichage des 30 premiers chiffres old <- par(no.readonly = TRUE) par(mfrow = c(5,6), mar = c(2,2,2,2)) for(i in 1:30) { y = as.matrix(train[i, 2:785]) dim(y) = c(28, 28) image(y[,ncol(y):1], axes = FALSE, col = gray(0:255/255)) } par(old) # chiffres moyens : superposition des images correspondant à un chiffre donné chiffres.moy <- matrix(NA,10,785) chiff.moy = function() { for (i in (0:9)) { chiffres.moy[i+1,1] <- i chiffres.moy[i+1,2:785] <- t(as.matrix(colMeans(train[which(train$label==i),2:785]))) } return(chiffres.moy) } chiffres.moy <- chiff.moy() # affichage des "chiffres moyens" old <- par(no.readonly = TRUE) par(mfrow = c(3,4), mar = c(2,2,2,2)) for(i in 1:10) { y = chiffres.moy[i, 2:785] dim(y) = c(28, 28) image(y[,ncol(y):1], axes = FALSE, col = gray(0:255/255)) } par(old) # score en fonction de la corrélation des pixels de l'image testée et de chaque chiffre moyen score = function(i) { dist.chiffr <- sapply(1:10, function(j) cor(as.numeric(test[i,-1]), chiffres.moy[j,2:785], method = "spearman")) return(which.max(dist.chiffr)-1) } system.time(pred.chif.moy <- Vectorize(score)(1:nrow(test))) # chiffre prédit par la corrélation maximale de l'observation aux 10 chiffres moyens pred.cm <- pred.chif.moy sum(pred.cm != test[,1]) / nrow(test) table(test[,1], pred.cm) # prédiction par maximisation de la corrélation entre les composantes principales des images et des composantes des 10 chiffres moyens X <- train[,-1] Y <- train[,1] Xreduced <- X/255 # ACP Xcov <- cov(Xreduced) pcaX <- prcomp(Xcov) # valeurs propres ev <- pcaX$sdev^2 # valeurs propres head(ev) head(ev/sum(ev)) head(cumsum(ev/sum(ev)), 200) # sommes cumulées des variances plot(pcaX) screeplot(pcaX) # les deux plots sont identiques # matrice d'apprentissage en entrée de la régression nbaxes <- 58 # matrice avec les 784 pixels en lignes les 58 axes en colonnes Xfinal <- as.matrix(Xreduced) %*% pcaX$rotation[,1:nbaxes] # données de test transformées par ACP testreduced <- test[,-1]/255 testreduced <- as.matrix(testreduced) %*% pcaX$rotation[,1:nbaxes] # images moyennes transformées par ACP chiffres.moy.reduced <- chiffres.moy[,-1] chiffres.moy.reduced <- data.frame(as.matrix(chiffres.moy.reduced) %*% pcaX$rotation[,1:nbaxes]) # score en fonction de la corrélation des pixels de l'image testée et de chaque chiffre moyen score = function(i) { dist.chiffr <- sapply(1:10, function(j) cor(as.numeric(testreduced[i,]), as.numeric(chiffres.moy.reduced[j,]), method = "pearson")) return(which.max(dist.chiffr)-1) } system.time(pred.chif.moy <- Vectorize(score)(1:nrow(test))) # chiffre prédit par la corrélation maximale de l'observation aux 10 chiffres moyens pred.cm <- pred.chif.moy sum(pred.cm != test[,1]) / nrow(test) table(test[,1], pred.cm) # Section 7.4 : Analyse discriminante linéaire ou quadratique # --------------------------------------------------------------------------------------------------------- # analyse discriminante linéaire library(MASS) system.time(lda <- lda(label ~ ., data=train)) # KO à cause des variables constantes dans les groupes varcst <- (1+c(1,2,3,4,5,6,7,8,9,10,11,12,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,53,54,55,56,57,58,83,84,85,86,112,113,141,142,169,477,561,645,646,672,673,674,700,701,702,728,729,730,731,755,756,757,758,759,760,781,782,783,784)) system.time(lda <- lda(label ~ ., data=train[,-varcst])) # prédiction par LDA et taux d'erreur en test pred.lda <- predict(lda, test)$class sum(pred.lda != test[,1]) / nrow(test) table(test[,1], pred.lda) # analyse discriminante quadratique library(MASS) # fonction jitter pour ajout de bruit aux données et suppression de la colinéarité train.J <- train train.J[,-1] <- apply(train.J[,-1], 2, jitter) system.time(qda <- qda(label ~ ., train.J)) pred.qda <- predict(qda, test)$class sum(pred.qda != test[,1]) / nrow(test) table(test[,1], pred.qda) # régression sur facteurs ACP # création ci-dessus des matrices Xfinal et testreduced trainreduced <- data.frame(Y, Xfinal) # remarque : la création d'un modèle directement à partir des vecteurs sans passer par un data frame pose un pb # dans la fonction predict, qui s'attend à un newdata de même longueur que l'échantillon d'apprentissage (60000 lignes) system.time(adq <- qda(Y ~ ., data=trainreduced)) # remarque : il faut que "testreduced" soit un data frame pour lui appliquer la fonction "predict" de MASS pred.qda <- predict(adq, newdata=data.frame(testreduced))$class # taux d'erreur en test sum(pred.qda != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.qda) # Section 7.5 : Régression logistique multinomiale # --------------------------------------------------------------------------------------------------------- library(h2o) localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, max_mem_size = '6g', ice_root="D:/Data/h2o", nthreads=4) # import de la base MNIST system.time(train <- h2o.importFile(path = "D:/Data/MNIST/mnist_train.csv", header=TRUE)) system.time(test <- h2o.importFile(path = "D:/Data/MNIST/mnist_test.csv", header=TRUE)) train$label <- as.factor(train$label) test$label <- as.factor(test$label) # régression multinomiale system.time(glmh2o <- h2o.glm(y = 1, x = 2:785, training_frame = train, family = "multinomial", alpha=0, lambda_search = TRUE, validation_frame = test)) system.time(glmh2o <- h2o.glm(y = 1, x = 2:785, training_frame = train, family = "multinomial", solver="L_BFGS", alpha=1, lambda = 0.0001, max_iterations = 500, validation_frame = test)) system.time(glmh2o <- h2o.glm(y = 1, x = 2:785, training_frame = train, family = "multinomial", solver="L_BFGS", alpha=0, lambda = 0.001, max_iterations = 500, validation_frame = test)) system.time(glmh2o <- h2o.glm(y = 1, x = 2:785, training_frame = train, family = "multinomial", solver="IRLSM", alpha=1, lambda = 0.0001, max_iterations = 50, validation_frame = test)) print(glmh2o) # application du modèle à l'échantillon de test prediction <- h2o.predict(glmh2o, newdata=test) pred.label <- max.col(prediction[,-1]) - 1 # ne pas prendre en compte la 1ère colonne qui contient l'étiquette sum(pred.label != as.vector(test$label)) / nrow(test) table(as.vector(test$label), pred.label) # matrice de confusion # stacking (code absent du livre) # nombre de validations croisées nfolds <- 5 # apprentissage et validation croisée d'une régression glmh2o <- h2o.glm(y = 1, x = 2:785, training_frame = train, family = "multinomial", alpha=0, lambda_search = TRUE, nfolds = nfolds, fold_assignment = "Modulo", keep_cross_validation_predictions = TRUE, seed = 235) # apprentissage et validation croisée d'un modèle de gradient boosting gbmh2o <- h2o.gbm(y = 1, x = 2:785, training_frame = train, ntrees = 300, max_depth = 10, min_rows = 1, learn_rate = 0.1, distribution= "multinomial", nfolds = nfolds, fold_assignment = "Modulo", keep_cross_validation_predictions = TRUE, seed = 235) # apprentissage et validation croisée d'une forêt aléatoire rfh2o <- h2o.randomForest(y = 1, x = 2:785, training_frame = train, ntrees = 300, max_depth = 30, mtries=20, nfolds = nfolds, fold_assignment = "Modulo", keep_cross_validation_predictions = TRUE, seed = 235) # apprentissage du modèle par stacking des modèles précédents ensemble <- h2o.stackedEnsemble(y = 1, x = 2:785, training_frame = train, model_id = "my_ensemble", base_models = list(glmh2o, gbmh2o, rfh2o)) h2o.confusionMatrix(ensemble, test) prediction <- h2o.predict(ensemble, newdata=test) pred.label <- max.col(prediction[,-1]) - 1 # ne pas prendre en compte la 1ère colonne qui contient l'étiquette sum(pred.label != as.vector(test$label)) / nrow(test) table(as.vector(test$label), pred.label) # matrice de confusion # modèle AutoML (code absent du livre) system.time(automl.model <- h2o.automl(y = 1, x = 2:785, training_frame = train, validation_frame = test, seed = 235, max_runtime_secs=1000, max_models=10)) summary(automl.model) automl.model@leaderboard h2o.confusionMatrix(automl.model, test) prediction <- h2o.predict(automl.model, newdata=test) pred.label <- max.col(prediction[,-1]) - 1 # ne pas prendre en compte la 1ère colonne qui contient l'étiquette sum(pred.label != as.vector(test$label)) / nrow(test) table(as.vector(test$label), pred.label) # matrice de confusion # Section 7.6 : Forêts aléatoires # --------------------------------------------------------------------------------------------------------- library(ranger) train[,1] <- as.factor(train[,1]) # la variable à expliquer doit être un facteur pour que ranger effectue un classement system.time(rg <- ranger(label ~ ., data=train, importance = "none", num.trees=100, mtry=30, replace=T, write.forest=T, seed=235)) rg # prédiction en apprentissage table(train[,1], rg$predictions) # prédiction en test pred.rg <- predict(rg, dat = test[,-1]) sum(pred.rg$predictions != test[,1]) / nrow(test) table(test[,1], pred.rg$predictions) # Section 7.7 : Extra-Trees # --------------------------------------------------------------------------------------------------------- #Sys.setenv(JAVA_HOME='C:/Program Files/Java/jdk1.8.0_181') options(java.parameters = "-Xmx3g") # 3g définit la mémoire allouée - 1 Gb n'est pas suffisant pour le MNIST library(extraTrees) cible <- as.factor(train[,1]) # la variable à expliquer doit être un facteur pour que les Extra-Trees effectuent un classement set.seed(235) system.time(et <- extraTrees(train[,-1], cible, ntree=300, mtry=10, numRandomCuts=3, nodesize=1, numThreads=1)) # prédiction en test pred.et <- predict(et, test[,-1], probability=F) sum(pred.et != test[,1]) / nrow(test) table(test[,1],pred.et) # Section 7.8 : Gradient boosting # --------------------------------------------------------------------------------------------------------- #library(devtools) #install_github('dmlc/xgboost', subdir='R-package') install.packages('xgboost') library(xgboost) # lecture des données train <- read.csv("D:/Data/MNIST/mnist_train.csv", header=TRUE) test <- read.csv("D:/Data/MNIST/mnist_test.csv", header=TRUE) # transformation data frames en matrices sparses library(Matrix) train.mx <- sparse.model.matrix(label~., train) test.mx <- sparse.model.matrix(label~., test) # transformation matrices sparses en matrices Xgb dtrain <- xgb.DMatrix(train.mx, label=train$label) dtest <- xgb.DMatrix(test.mx, label=test$label) # modèle avec double boucle sur paramètres êta et colsample_bytree pour optimisation i <- seq(0.05,0.25,by=0.05) # 5 valeurs j <- seq(0.2,0.5,by=0.05) # 7 valeurs f <- function(i,j) { set.seed(123) #train.gdbt <- xgb.train(params=list(booster = "gbtree", objective="multi:softmax", num_class=10, eval_metric="mlogloss", max_depth=5, subsample=1, eta=i, colsample_bytree=j), # data=dtrain, nrounds=100, nthread = 4, verbose=0, watchlist=list(eval=dtest, train=dtrain)) train.gdbt <- xgb.train(params=list(booster = "gbtree", objective="multi:softmax", num_class=10, eval_metric="mlogloss", max_depth=5, subsample=1, eta=i, colsample_bytree=j), data=dtrain, nrounds=1000, early_stopping_rounds =10, nthread = 4, verbose=0, watchlist=list(eval=dtest)) pred.gbm <- predict(train.gdbt, newdata=dtest) return(sum(pred.gbm != test$label) / nrow(test)) } system.time(k <- outer(i,j,Vectorize(f))) k which(k==min(k),arr.ind=TRUE) (eta <- i[which(k==min(k),arr.ind=TRUE)[1]]) (colsample_bytree <- j[which(k==min(k),arr.ind=TRUE)[2]]) # modélisation set.seed(123) system.time(train.gdbt <- xgb.train(params=list(booster = "gbtree", objective="multi:softmax", num_class=10, eval_metric="mlogloss", eta=0.15, max_depth=5, subsample=1, colsample_bytree=0.2), data=dtrain, nrounds=1000, nthread = 4, verbose=0)) system.time(train.gdbt <- xgb.train(params=list(booster = "gbtree", objective="multi:softmax", num_class=10, eval_metric="mlogloss", eta=0.15, max_depth=5, subsample=1, colsample_bytree=0.2), data=dtrain, nrounds=1000, nthread = 4, watchlist=list(eval=dtest), early_stopping_rounds =10, verbose=1)) pred.gbm <- predict(train.gdbt, newdata=dtest) sum(pred.gbm != test$label) / nrow(test) table(test$label, pred.gbm) # variable importance plot mat <- xgb.importance(feature_names = colnames(train), model = train.gdbt) xgb.plot.importance(importance_matrix = mat[1:20]) # 20 variables les plus discriminantes # boosting avec h2o system.time(gbmh2o <- h2o.gbm(x = 2:785, y = 1, training_frame = train, ntrees = 300, max_depth = 10, min_rows = 1, learn_rate = 0.1, distribution= "multinomial", validation_frame = test)) print(gbmh2o) # application du modèle à l'échantillon de test prediction <- h2o.predict(gbmh2o, newdata=test) pred.label <- max.col(prediction[,-1]) - 1 # ne pas prendre en compte la 1ère colonne qui contient l'étiquette sum(pred.label != as.vector(test$label)) / nrow(test) table(as.vector(test$label), pred.label) # matrice de confusion # Section 7.9 : Machines à vecteurs de support (Support Vector Machines) # --------------------------------------------------------------------------------------------------------- # SVM sur axes factoriels (calculés précédemment en section 7.3) library(e1071) # SVM avec noyau polynomial svmpol <- tune.svm(Xfinal, factor(Y), kernel="polynomial", degree=3, gamma=10^(-3:0), cost=10^(-1:2), coef0=seq(0,1,by=0.1), validation.x = testreduced, validation.y = factor(test[,1]), tunecontrol = tune.control(sampling = "fix", fix=1)) svmpol$best.model summary(svmpol) system.time(svmpol <- svm(Xfinal, factor(Y), kernel="polynomial", degree=2, coef0=0.1, gamma=0.01, cost=10, probability=TRUE)) system.time(svmpol <- svm(Xfinal, factor(Y), kernel="polynomial", degree=3, coef0=0.2, gamma=0.01, cost=10, probability=TRUE)) print(svmpol) # prédiction pred.svm <- predict(svmpol, newdata=testreduced, decision.values=T) # taux d'erreur en test sum(pred.svm != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.svm) # SVM avec noyau radial system.time(svmrad <- svm(Xfinal, factor(Y), kernel="radial", probability=TRUE, gamma=0.01, cost=10)) print(svmrad) # prédiction pred.svm <- predict(svmrad, newdata=testreduced, decision.values=T) # taux d'erreur en test sum(pred.svm != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.svm) # Section 7.10 : Perceptron à une couche cachée # --------------------------------------------------------------------------------------------------------- # perceptron : modélisation directe (sans passer par une ACP) # avec une variable à expliquer facteur : pas besoin de préciser une sortie softmax library(nnet) set.seed(235) system.time(rn <- nnet(factor(label) ~ ., data=train, size=1, maxit=100)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=1, decay=1, maxit=100)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=2, maxit=100, MaxNWts =10000)) system.time(rn <- nnet(factor(label) ~ ., data=train, size=4, maxit=100, MaxNWts =10000, decay=0.5)) (nb_poids <- (784+10)*1 + (10+1)) # calcul du nb de pondérations du réseau # prédiction pred.rn <- predict(rn, newdata=test[,-1], type="class") head(pred.rn) # taux d'erreur en test sum(pred.rn != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.rn) # perceptron sur facteurs ACP avec sortie softmax set.seed(235) trainDataset <- train[sample(nrow(train)),] # mélange des images X <- trainDataset[,-1] Y <- trainDataset[,1] Xreduced <- X/255 # ACP Xcov <- cov(Xreduced) pcaX <- prcomp(Xcov) # matrice d'apprentissage en entrée du réseau de neurones nbaxes <- 58 Xfinal <- as.matrix(Xreduced) %*% pcaX$rotation[,1:nbaxes] Y <- class.ind(Y) # données de test transformées par ACP testreduced <- test[,-1]/255 testreduced <- as.matrix(testreduced) %*% pcaX$rotation[,1:nbaxes] # PMP sur axes factoriels set.seed(235) system.time(rn <- nnet(Xfinal, Y, size=23, softmax=TRUE, maxit=100, MaxNWts =10000)) system.time(rn <- nnet(Xfinal, Y, size=23, softmax=TRUE, maxit=100, MaxNWts =10000, decay=1)) system.time(rn <- nnet(Xfinal, Y, size=100, softmax=TRUE, maxit=300, MaxNWts =30000, decay=1)) rn # prédiction pred.rn <- predict(rn, newdata=testreduced, type="class") # taux d'erreur en test sum(pred.rn != test[,1]) / nrow(test) # matrice de confusion table(test[,1], pred.rn) # Section 7.11 : Réseau de neurones « deep learning » de H2O # --------------------------------------------------------------------------------------------------------- library(h2o) localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, max_mem_size = '6g', ice_root="D:/Data/h2o", nthreads=4) # import de la base MNIST system.time(train <- h2o.importFile(path = "D:/Data/MNIST/mnist_train.csv", header=TRUE)) system.time(test <- h2o.importFile(path = "D:/Data/MNIST/mnist_test.csv", header=TRUE)) train$label <- as.factor(train$label) test$label <- as.factor(test$label) # modèles de réseau de neurones dlh2o <- h2o.deeplearning(x = 2:785, y = 1, training_frame = train, validation_frame = test, activation = "Tanh", input_dropout_ratio = 0.2, hidden = c(50,50,50), epochs = 10) dlh2o <- h2o.deeplearning(x = 2:785, y = 1, training_frame = train, validation_frame = test, activation = "Rectifier", input_dropout_ratio = 0.2, hidden = c(50,50,50), epochs = 10) dlh2o <- h2o.deeplearning(x = 2:785, y = 1, training_frame = train, validation_frame = test, activation = "Rectifier", max_w2 = 10, l2 = 1e-3, input_dropout_ratio = 0.2, hidden = c(50,50,50), epochs = 10) dlh2o <- h2o.deeplearning(x = 2:785, y = 1, training_frame = train, validation_frame = test, activation = "Rectifier", max_w2 = 20, l2 = 1e-3, input_dropout_ratio = 0.2, hidden = c(500,500,500), epochs = 20) system.time(dlh2o <- h2o.deeplearning(x = 2:785, y = 1, training_frame = train, validation_frame = test, activation = "Rectifier", l1 = 1e-5, max_w2 = 5, input_dropout_ratio = 0.2, hidden = c(1024,1024,2048), epochs = 20)) # autre bon modèle system.time(dlh2o <- h2o.deeplearning(x = 2:785, y = 1, training_frame = train, validation_frame = test, activation = "RectifierWithDropout", max_w2 = 20, l2 = 1e-3, input_dropout_ratio = 0.2, hidden_dropout_ratios = c(0.5,0.5,0.5), hidden = c(500,500,500), balance_classes = FALSE, seed=123, epochs = 20)) print(dlh2o) # application du modèle à l'échantillon de test prediction <- h2o.predict(dlh2o, newdata=test) pred.label <- max.col(prediction[,-1]) - 1 # ne pas prendre en compte la 1ère colonne qui contient l'étiquette sum(pred.label != as.vector(test$label)) / nrow(test) table(as.vector(test$label), pred.label) # matrice de confusion # arrêt du cluster h2o.shutdown(prompt=FALSE) # Section 7.13 : Réseau de neurones convolutif # --------------------------------------------------------------------------------------------------------- library(mxnet) #detach(package:h2o) # fichiers MNIST train <- read.csv("D:/Data/MNIST/mnist_train.csv", header=TRUE) test <- read.csv("D:/Data/MNIST/mnist_test.csv", header=TRUE) # préparation des données pour mxnet train.x <- train[,-1] train.y <- train[,1] train.x <- t(train.x/255) test.x <- test[,-1] test.x <- t(test.x/255) # 784 lignes et 10000 colonnes, valeurs entre 0 et 1 # transformation des matrices en arrays (nécessaire pour l'opérateur de convolution de mxnet) : train.array <- train.x dim(train.array) <- c(28, 28, 1, ncol(train.x)) test.array <- test.x dim(test.array) <- c(28, 28, 1, ncol(test.x)) # 2 couches de convolution avec pooling recouvrant et 1 dropout # entrée data <- mx.symbol.Variable('data') # 1ère couche de convolution conv1 <- mx.symbol.Convolution(data=data, kernel=c(5,5), num_filter=20, workspace=4096) tanh1 <- mx.symbol.Activation(data=conv1, act_type="tanh") pool1 <- mx.symbol.Pooling(data=tanh1, pool_type="max", kernel=c(3,3), stride=c(2,2)) # 2e couche de convolution conv2 <- mx.symbol.Convolution(data=pool1, kernel=c(5,5), num_filter=60, workspace=4096) tanh2 <- mx.symbol.Activation(data=conv2, act_type="tanh") pool2 <- mx.symbol.Pooling(data=tanh2, pool_type="max", kernel=c(3,3), stride=c(2,2)) # couche complètement connectée flatten <- mx.symbol.Flatten(data=pool2) drp <- mx.symbol.Dropout(data=flatten, p=0.5) fc1 <- mx.symbol.FullyConnected(data=drp, num_hidden=600) tanh3 <- mx.symbol.Activation(data=fc1, act_type="tanh") # couche de sortie fc2 <- mx.symbol.FullyConnected(data=tanh3, num_hidden=10) # loss lenet <- mx.symbol.SoftmaxOutput(data=fc2) # apprentissage device <- mx.cpu() mx.set.seed(123) system.time(model <- mx.model.FeedForward.create(lenet, X=train.array, y=train.y, ctx=device, num.round=30, array.batch.size=100, learning.rate=0.1, momentum=0.65, wd=0.00001, eval.metric=mx.metric.accuracy, epoch.end.callback=mx.callback.log.train.metric(100))) # application du modèle preds <- predict(model, test.array) pred.label <- max.col(t(preds)) - 1 (nrow(test)-sum(test$label == pred.label))/nrow(test) # taux d'erreur table(test$label,pred.label) # définition d'une fonction "mode" Mode <- function(x) { if (is.numeric(x)) { x_table <- table(x) return(as.numeric(names(x_table)[which.max(x_table)])) } } # comité de réseaux à taux d'apprentissage et moment variant aléatoirement comitcnn <- function(nc=10, nr=10, lr=0.1, mm=0.7) { pred.list <- list() for (i in 1:nc) { model <- mx.model.FeedForward.create(lenet, X=train.array, y=train.y, ctx=device, num.round=nr, array.batch.size=128, initializer = mx.init.uniform(0.05), learning.rate=lr + runif(n=1, min=-0.05, max=0.05), momentum=mm + runif(n=1,min=-0.1,max=0.1), wd=0.00001, eval.metric=mx.metric.accuracy, epoch.end.callback=mx.callback.log.train.metric(100)) preds <- predict(model, test.array) pred.list[[i]] <- max.col(t(preds)) - 1 } cnn.res <- do.call("cbind", pred.list) pred.label <- apply(cnn.res,1,Mode) print(table(test$label,pred.label)) erreur <- (nrow(test)-sum(diag(table(test$label,pred.label))))/nrow(test) print(erreur) return(pred.label) } system.time(comitcnn(10,10,0.1,0.7)) # --------------------------------------------------------------------------------------------------------- # Chapitre 8 : Big Data avec R # --------------------------------------------------------------------------------------------------------- # Section 8.2 : Langage interprété vs compilé # --------------------------------------------------------------------------------------------------------- library(microbenchmark) microbenchmark(calcul1={ I = 0 while (I < 100000) { I = I + 1 } }, calcul2={ I = 0 while (I < 100000) { ((((((I = I + 1)))))) } }, times=100, unit="ms") # le passage par une fonction gomme les écarts dus aux parenthèses car cette fonction est compilée library(compiler) f <- function(n){ I = 0 while (I < n) { I = I + 1 } return(I) } g <- function(n){ I = 0 while (I < n) { ((((((I = I + 1)))))) } return(I) } f1 <- cmpfun(f) g1 <- cmpfun(g) microbenchmark(f(100000), g(100000), f1(100000), g1(100000), times=100, unit="ms") # Section 8.3 : Gestion de la mémoire par R # --------------------------------------------------------------------------------------------------------- # allocation RAM sous Windows (pas sous Linux) memory.size(max=NA) memory.limit(size=6096) # taille des objets rev(sort(sapply(ls(),function(x){object.size(get(x))}))) (RAM <- rev(sort(sapply(ls(),function(x){object.size(get(x))})))) barplot(RAM) sum(gc()[,1]) # ménage de la mémoire rm(list=ls()) # pour tout effacer do_not_remove <- c("object1", "object2") rm(list = ls()[!(ls() %in% do_not_remove)]) gc() gc(reset=T) objet <- seq(100000) gc() rm(objet) gc(reset=T) objet <- as.list(seq(100000)) gc() # effet de garbage collection sur allocation mémoire (voir le gestionnaire des tâches Windows) x <- rep (0, 2^27) object.size(x) rm(x) gc() # Section 8.3 : Profiling de la mémoire vive # --------------------------------------------------------------------------------------------------------- # profiling de la mémoire avec fonction Rprof Rprof("profile1.out") example(glm) Rprof(NULL) summaryRprof("profile1.out") summaryRprof("profile1.out")$by.total # diagnostic ligne par ligne Rprof("profile2.out", line.profiling=TRUE) example(glm) Rprof(NULL) summaryRprof("profile2.out", lines = "show") profil <- summaryRprof("profile2.out", lines = "show")$by.total (result <- cbind("function"=rownames(profil), profil)) # profiling de la mémoire avec package profr #install.packages('profr') library(profr) prof <- profr(example(glm)) prof plot(prof) head(prof[rev(order(prof$time)),],20) # profiling de la mémoire avec package lineprof # le package lineprof est absent du CRAN library(devtools) install_github("hadley/lineprof") library(lineprof) # N.B. L'option "torture = TRUE" multiplie le temps de calcul par 10 à 100 prof <- lineprof(example(glm), torture = FALSE) prof shine(prof) # Section 8.3 : SOAR # --------------------------------------------------------------------------------------------------------- install.packages('SOAR') r <- data.frame(a=rnorm(10,2,.5),b=rnorm(10,3,.5)) r ls() library(SOAR) Store(r) ls() mean(r[,1]) ls() r$c = rnorm(10,4,.5) ls() r head(r,2) # Section 8.4 : Limites de la mémoire # --------------------------------------------------------------------------------------------------------- x <- integer(2^31-1) x <- numeric(2^31-1) x <- matrix(0,2^31,1) x <- integer(2^53) x <- matrix(0,2^30,2) x <- matrix(0,2^30,2^22) x <- matrix(0,2^30,2^23) library(ff) x <- ff(0,length = 2^31) x <- ff(0,length = 2^31-1) x <- ff(0,dim=c(2^31-1,1)) dim(x) x <- ff(0,dim=c(2^31,1)) x <- ff(0,dim=c(2^30,2)) # Section 8.5.1 : Programmation dans R : taille des objets influencées par type # --------------------------------------------------------------------------------------------------------- object.size(matrix(0,100000,30)) object.size(matrix(as.integer(0),100000,30)) library('Matrix') object.size(Matrix(0,100000,30)) # Section 8.5.2 : Les boucles dans R # --------------------------------------------------------------------------------------------------------- # boucle avec allocation préalable time <- proc.time() x <- seq(1:1000000) res <- numeric(1000000) for (i in 1:1000000) res[i] = exp(x[i]) proc.time() - time # boucle sans allocation préalable time <- proc.time() x <- seq(1:1000000) res <- numeric() # équivalent à numeric(0) : crée une variable numérique mais sans donnée for (i in 1:1000000) res[i] = exp(x[i]) proc.time() - time # effet de l'allocation d'espace dans la création d'une matrice system.time({set.seed(123) ; ans1 = numeric() ; for(i in 1:10000) ans1 = cbind(ans1, rnorm(10))}) system.time({set.seed(123) ; ans2 = matrix(NA, 10, 10000) ; for(i in 1:10000) ans2[,i] = rnorm(10)}) system.time({set.seed(123) ; ans3 <- sapply(1:10000, function(i) rnorm(10))}) system.time({set.seed(123) ; ans1 = numeric() ; for(i in 1:100000) ans1 = cbind(ans1, rnorm(10))}) system.time({set.seed(123) ; ans2 = matrix(NA, 10, 100000) ; for(i in 1:100000) ans2[,i] = rnorm(10)}) system.time({set.seed(123) ; ans3 <- sapply(1:100000, function(i) rnorm(10))}) dim(ans1) identical(ans1,ans2, ans3) # outer n <- 10000 i <- rnorm(n) j <- rnorm(n) time <- proc.time() z <- matrix(0,n,n) for (x in 1:n) { for (y in 1:n) { z[x,y] <- i[x]*j[y] } } proc.time() - time system.time(z1 <- outer(i,j,"*")) dim(z1) identical(z,z1) # Section 8.5.3 : Vectorisation # --------------------------------------------------------------------------------------------------------- lst <- list(a=c(1,2,4), b=c(1,4,16)) lst lapply(lst,sqrt) lapply(lst,mean) sapply(lst,sqrt) sapply(lst,mean) system.time(res <- lapply(1:1000000, exp)) system.time(res <- vapply(1:1000000, exp, FUN.VALUE=0)) i39 <- lapply(3:9, seq) sapply(i39, fivenum) vapply(i39, fivenum,c(Min. = 0, "1st Qu." = 0, Median = 0, "3rd Qu." = 0, Max. = 0)) # une régression linéaire par groupe : avec lapply data(iris) summary(lm(Sepal.Length ~ Petal.Length, iris)) system.time(lst.model <- lapply(split(iris, iris$Species), function(d) lm(Sepal.Length ~ Petal.Length, d))) system.time(lst.model <- lapply(split(iris, iris$Species), lm, formula=Sepal.Length ~ Petal.Length)) lst.model lapply(lst.model, coef) # extraction des coefficients g <- function(m) coef(m)[1] h <- function(m) coef(m)[2] (coeff <- data.frame(cste=unlist(lapply(lst.model, g)), pente=unlist(lapply(lst.model, h)), row.names = names(lst.model))) # prédictions mapply('predict', lst.model, split(iris, iris$Species)) # application d'un modèle à une liste de data frames training <- data.frame(y=rnorm(10), x1=rnorm(10), x2=rnorm(10)) model <- lm(y~., data=training) list.df <- list() for (i in 1:5){ list.df[[i]] <- data.frame(x1=rnorm(10), x2=rnorm(10)) } predictions <- lapply(list.df, function(x) {predict(model, newdata=x)}) # Section 8.5.4 : Appel de code C++ dans R # --------------------------------------------------------------------------------------------------------- # on doit utiliser le programme Rtools pour compiler un package binaire sous Windows # à noter qu'à partir de la version 4.0.0 de R, Rtools est remplacé par Rtools40 (https://cran.r-project.org/bin/windows/Rtools/) # variables d'environnement Windows à créer pour indiquer à R où est installé Rtools Sys.setenv(PATH = paste("D:/Rtools/bin", Sys.getenv("PATH"), sep=";")) Sys.setenv(BINPREF = "D:/Rtools/mingw_$(WIN)/bin/") library(Rcpp) evalCpp("1+1") # exemple 1 cppFunction('NumericVector expC(NumericVector x) { int n = x.size(); NumericVector out(n); for(int i = 0; i < n; ++i) { out[i] = exp(x[i]); } return out; }') expC(1:100) # exemple 2 cppFunction('NumericVector expC2(NumericVector x) { return exp(x); }') expC2(1:100) # Section 8.5.5 : Comparaison des temps de calcul # --------------------------------------------------------------------------------------------------------- nb <- 1000000 x <- seq(1:nb) res <- numeric(nb) f <- function(n) { exp(n) } library(compiler) g <- cmpfun(f) library(snow) clus <- makeCluster(4) library(parallel) # pour fonction mclapply library(microbenchmark) compare <- microbenchmark(exp(1:nb), Vectorize(exp)(1:nb), sapply(1:nb, function(x) exp(x)), sapply(1:nb, f), sapply(1:nb, g), sapply(1:nb, exp), lapply(1:nb, exp), vapply(1:nb, exp, FUN.VALUE=0), parSapply(clus,1:nb,exp), parLapply(clus,1:nb,exp), mclapply(1:nb,exp), parSapply(clus,seq(1,nb,100),function(n) exp(n+(0:99))), parLapply(clus,seq(1,nb,100),function(n) exp(n+(0:99))), mclapply(seq(1,nb,100),function(n) exp(n+(0:99))), for (i in 1:nb) res[i] = exp(x[i]), expC(1:nb),expC2(1:nb), source("sourceR_notepad.R"), times=100, unit="ms") compare library(ggplot2) autoplot(compare) # Section 8.5.6 : Méthodes de Monte-Carlo # --------------------------------------------------------------------------------------------------------- # calcul de PI par intégrale de Monte-Carlo N <- 10000000 f = function(N){ hits = 0 for(i in 1:N) { u1 = runif(1,-1,1); u2 = runif(1,-1,1) if((u1^2 + u2^2) < 1) hits = hits + 1 } return(4*hits/N) } system.time(PI <- f(N)) PI (PI-pi)/pi g = function(N){ hits = sum(runif(N,-1,1)^2 < (1- runif(N,-1,1)^2)) return(4*hits/N) } system.time(PI1 <- g(N)) PI1 (PI1-pi)/pi # Section 8.5.7 : Calcul rapide de dates # --------------------------------------------------------------------------------------------------------- library(data.table) system.time(Test <- fread("D:/Data/vols/2008.csv")) library(stringr) system.time(Test[,':='(Date = paste0(Year,str_pad(Month,2,pad = "0"),str_pad(DayofMonth,2,pad = "0")), Heure = paste0(substr(str_pad(CRSDepTime,4,pad = "0"),1,2),substr(str_pad(CRSDepTime,4,pad = "0"),3,4),"00"))]) head(Test[,.(Date,Heure)]) system.time(Test[,':='(Temps = as.POSIXct(paste(Date,Heure), format='%Y%m%d %H%M%S'))]) head(Test[,.(Date,Heure,Temps)]) Test[,Temps:=NULL] # suppression de la variable Temps library(fasttime) system.time(Test[,':='(Temps = fastPOSIXct(paste(paste0(substr(Date,1,4),"-", substr(Date,5,6),"-", substr(Date,7,8)), paste0(substr(Heure,1,2),":", substr(Heure,3,4),":", substr(Heure,5,6))), "GMT"))]) head(Test[,.(Date,Heure,Temps)]) # Section 8.5.8 : Analyse factorielle # --------------------------------------------------------------------------------------------------------- # packages benchmarkés library(FactoMineR) library(MASS) library(gmodels) library(ade4) library(microbenchmark) library(ggplot2) # jeu de données avec des variables quantitatives n <- 10000 # nombre de lignes p <- 100 # nombre de colonnes myVars <- paste0("V",1:p) X <- data.frame(matrix(ncol = length(myVars), nrow = n)) names(X) <- myVars for (i in myVars) { X[,i] <- rnorm(n) } dim(X) # benchmark des ACP compare <- microbenchmark(Facto=PCA(X, ncp = 99, graph = F), prcmp=prcomp(X, scale. = TRUE, rank = 99), princmp=princomp(X, cor = TRUE), fastpca=fast.prcomp(X, scale = TRUE), ade4=dudi.pca(X, scannf = F, nf = 99), times=30, unit="ms") compare autoplot(compare) # jeu de données avec des variables qualitatives n <- 10000 # nombre de lignes p <- 100 # nombre de colonnes myVars <- paste0("V",1:p) Y <- data.frame(matrix(ncol = length(myVars), nrow = n)) names(Y) <- myVars for (i in myVars) { Y[,i] <- factor(LETTERS[(sample(n,n)%% 26)+1]) } dim(Y) # benchmark des ACM compare <- microbenchmark(Facto=MCA(Y, ncp = 99, graph = F), MASS=mca(Y, nf = 99), ade4=dudi.acm(Y, scannf = F, nf = 99), times=30, unit="ms") compare autoplot(compare) # Section 8.6.1 : Parallélisation des traitements # --------------------------------------------------------------------------------------------------------- library(parallel) detectCores() detectCores(logical=FALSE) # uniquement les coeurs physiques library(snow) clus <- makeCluster(4) time <- proc.time() res <- parSapply(clus,1:1000000,exp) proc.time() - time head(res) time <- proc.time() res <- parLapply(clus,1:1000000,exp) proc.time() - time head(res) head(unlist(res)) # réduction du nombre de tâches distribuées et donc du temps de calcul library(microbenchmark) (compare <- microbenchmark(T10000=parLapply(clus, seq(1,1000000,100), function(n) exp(n+(0:99))), T1000=parLapply(clus, seq(1,1000000,1000), function(n) exp(n+(0:999))), T100=parLapply(clus, seq(1,1000000,10000), function(n) exp(n+(0:9999))), T10=parLapply(clus, seq(1,1000000,100000), function(n) exp(n+(0:99999))), T1=parLapply(clus, seq(1,1000000,1000000), function(n) exp(n+(0:999999))), times=100, unit="ms")) # passage d'un paramètre res <- parLapply(clus,1:1000000,function(x) exp(2*x)) head(unlist(res)) k <- 2 clusterExport(clus,"k") res <- parLapply(clus,1:1000000,function(x) exp(k*x)) head(unlist(res)) # foreach library(doSNOW) cl <- makeCluster(4) registerDoSNOW(clus) library(foreach) foreach (i=1:3) %do% { sqrt(i) } foreach (i=1:3, .combine=c) %do% { sqrt(i) } foreach (i=1:3, .combine=cbind) %do% { sqrt(i) } foreach (i=1:3, .combine=rbind) %do% { sqrt(i) } foreach (i=1:3, .combine='+') %do% { sqrt(i) } # forêts aléatoires parallélisées library(randomForest) # noms des variables myVariableNames <- c("Comptes","Duree_credit","Historique_credit","Objet_credit", "Montant_credit","Epargne","Anciennete_emploi","Taux_effort","Situation_familiale", "Garanties","Anciennete_domicile","Biens","Age","Autres_credits","Statut_domicile", "Nb_credits","Type_emploi","Nb_pers_charge","Telephone","Etranger","Cible") # jeu de données credit = read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data", h=FALSE, col.names=myVariableNames) credit$Cible <- as.factor(credit$Cible) # version sans foreach set.seed(235) # pour reproductibilité system.time(rf <- randomForest(Cible ~ ., data=credit, importance=TRUE, ntree=10000, mtry=3, replace=T, keep.forest=T, nodesize=5)) # appel foreach sans parallélisation library(foreach) ptm <- proc.time() rf <- foreach(ntree=rep(2500, 4), .combine=combine) %do% randomForest(Cible ~ ., data=credit, importance=TRUE, ntree=ntree, mtry=3, replace=T, keep.forest=T, nodesize=5) proc.time() - ptm # appel foreach avec parallélisation library(doSNOW) clus <- makeCluster(4) registerDoSNOW(clus) library(parallel) clusterSetRNGStream(clus,iseed=235) # pour reproductibilité ptm <- proc.time() rf <- foreach(ntree=rep(2500, 4), .combine=combine, .packages='randomForest') %dopar% randomForest(Cible ~ ., data=credit, importance=TRUE, ntree=ntree, mtry=3, replace=T, keep.forest=T, nodesize=5) proc.time() - ptm # fin parallélisation stopCluster(clus) # Section 8.6.2 : Parallélisation et génération de nombres aléatoires # --------------------------------------------------------------------------------------------------------- # C’est le générateur de Mersenne-Twister (dans la version de Matsumoto and Nishimura de 1998, et non la version améliorée de 2002) qui est utilisé par défaut dans R : RNGkind() # On crée une graine pour la reproductibilité des tirages aléatoires : seed <- 123 # On crée une fonction qui renvoie un nombre tiré selon la loi uniforme, exécutée nsim fois : nsim <- 100 testRNG <- function(i) { res <- runif(1) return(res) } # On voit la nécessité de spécifier une graine avant chaque tirage : set.seed(seed) res1 <- sapply(1:nsim, testRNG) res2 <- sapply(1:nsim, testRNG) identical(res1,res2) set.seed(seed) res1 <- sapply(1:nsim, testRNG) set.seed(seed) res2 <- sapply(1:nsim, testRNG) identical(res1,res2) # De façon plus simple : set.seed(seed) res1 <- runif(nsim) set.seed(seed) res2 <- runif(nsim) identical(res1,res2) # Nous voulons à présent paralléliser ce traitement et nous créons un cluster (le nombre de cœurs est automatiquement détecté) : library(parallel) clus <- makeCluster(detectCores()) # Il faut spécifier une graine et pouvoir générer un flux de nombres aléatoires dans les différents cœurs du cluster, mais la commande précédente ne le permet pas : set.seed(seed) res1 <- parSapply(clus, 1:nsim, testRNG) set.seed(seed) res2 <- parSapply(clus, 1:nsim, testRNG) identical(res1,res2) # Il faut utiliser la fonction suivante, qui permet de générer un flux de nombres aléatoires dans l’ensemble des cœurs du cluster, # en évitant de voir les mêmes séquences dans plusieurs cœurs, mais en évitant aussi d’être dépendant de la structure du cluster : clusterSetRNGStream(clus, seed) res1 <- parSapply(clus, 1:nsim, testRNG) clusterSetRNGStream(clus, seed) res2 <- parSapply(clus, 1:nsim, testRNG) identical(res1,res2) # Si l’on veut retrouver la même séquence avec la fonction sapply précédente, il faut spécifier le générateur de L'Ecuyer (1999), # qui est celui utilisé par la fonction clusterSetRNGStream et non le générateur de Mersenne-Twister utilisé dans R : RNGkind("L'Ecuyer-CMRG") set.seed(seed) res3 <- sapply(1:nsim, testRNG) clus <- makeCluster(1) clusterSetRNGStream(clus, seed) res1b <- parSapply(clus, 1:nsim, testRNG) identical(res1b,res3) # On peut revenir au générateur par défaut : library(randtoolbox) set.generator("default") RNGkind() # Section 8.8.2 : Le package data.table # --------------------------------------------------------------------------------------------------------- system.time(v2008 <- read.csv(file="D:/Data/vols/2008.csv", header=TRUE)) system.time(v2008 <- read.csv(file="D:/Data/vols/2008.csv", header=TRUE, stringsAsFactors=FALSE)) time <- proc.time() debut <- read.table(file="D:/Data/vols/2008.csv", header=TRUE, sep=',', nrows = 1000) classes <- sapply(debut, class) v2008 <- read.table(file="D:/Data/vols/2008.csv", header=TRUE, sep=',', colClasses = classes) proc.time() - time library(data.table) getDTthreads(verbose=TRUE) setDTthreads(threads = 0) getDTthreads() system.time(v2008.dt <- fread("D:/Data/vols/2008.csv")) all.equal(v2008, v2008.dt, check.attributes = F) # changement de nom d'une variable setnames(v2008.dt,"UniqueCarrier","Carrier") # suppression de variables v2008.dt[, c("Cancelled","CancellationCode") := NULL] ## création d’une nouvelle variable v2008.dt[, Age := v2008.dt[,Year]*12 + v2008.dt[,Month] - … ] # calcul du nombre de lignes de la data table v2008.dt[,.N] # dernière ligne de la table v2008.dt[.N,] # 1ère et dernière ligne pour chaque groupe défini par UniqueCarrier v2008.dt[,.SD[c(1,.N)], by=UniqueCarrier] # application d'une fonction pour chaque groupe défini par UniqueCarrier v2008.dt[,lapply(.SD, summary), by=UniqueCarrier] # sélection de lignes system.time(ans2 <- v2008[v2008$UniqueCarrier=="AA" & v2008$Origin=="JFK",]) # vector scan system.time(ans1 <- v2008.dt[UniqueCarrier=="AA" & Origin=="JFK",]) # mauvaise utilisation du data.table (avec vector scan) system.time(setkey(v2008.dt,UniqueCarrier,Origin)) # création d'un index system.time(ans1 <- v2008.dt[UniqueCarrier=="AA" & Origin=="JFK",]) # vector scan semble légèrement plus rapide avec index system.time(ans1 <- v2008.dt[list("AA","JFK")]) # recherche binaire system.time(ans1 <- v2008.dt[.("AA","JFK")]) # recherche binaire identical(ans1$FlightNum, ans2$FlightNum) # comparaison des temps d'agrégation par data.table, tapply et summarise system.time(ss <- v2008.dt[,mean(Distance),by=TailNum]) setkey(v2008.dt,TailNum) # création d'un index system.time(ss <- v2008.dt[,mean(Distance),by=TailNum]) system.time(tt <- tapply(v2008$Distance,v2008$TailNum,mean)) head(ss) head(tt) identical(as.vector(tt), ss$V1) # FALSE car l'ordre des éléments agrégés n'est pas le même identical(sort(as.vector(tt)),sort(ss$V1)) # TRUE # avec package dplyr library(dplyr) system.time(uu <- summarise(group_by(v2008.dt,TailNum),mean(Distance))) head(uu) all.equal(ss, as.data.table(uu), check.attributes = F) # agrégation comme ci-dessus avec ajout du nombre de lignes system.time(ss <- v2008.dt[,.(.N, distance.moy=mean(Distance)), by=TailNum]) head(ss) # agrégation sur plusieurs critères system.time(ss <- v2008.dt[,.(.N, distance.moy=mean(Distance)), by=.(UniqueCarrier,Origin)]) system.time(setkey(v2008.dt,UniqueCarrier,Origin)) # création d'un index system.time(ss <- v2008.dt[,.(.N, distance.moy=mean(Distance)), by=.(UniqueCarrier,Origin)]) head(ss) # avec dplyr library(dplyr) system.time(uu <- summarise(group_by(v2008.dt,UniqueCarrier,Origin), N=n(), distance.moy=mean(Distance))) head(uu) all.equal(ss, as.data.table(uu), check.attributes = F) # avec multidplyr #devtools::install_github("hadley/multidplyr") library(multidplyr) cluster <- new_cluster(8) system.time(vv <- v2008.dt %>% group_by(UniqueCarrier,Origin) %>% partition(cluster) %>% summarize(N=n(), distance.moy=mean(Distance)) %>% collect()) head(vv) all.equal(ss, as.data.table(vv), check.attributes = F) # KO all.equal(ss, as.data.table(vv)[order(vv$UniqueCarrier, vv$Origin),], check.attributes = F) # OK # distance moyenne d'un vol v2008.dt[,.(distance.moy=mean(Distance), retard.min=min(ArrDelay,na.rm=T), retard.max=max(ArrDelay,na.rm=T))] # fusion de data tables (dt1 <- data.table(A = letters[1:10], X = 1:10, key = "A")) (dt2 <- data.table(A = letters[5:14], Y = 1:10, key = "A")) merge(dt1, dt2) dt1[dt2, nomatch=0] merge(dt1, dt2, all = TRUE) merge(dt1, dt2, all.x = TRUE) dt2[dt1] merge(dt1, dt2, all.y = TRUE) dt1[dt2] library(microbenchmark) microbenchmark(merge(dt1, dt2), dt1[dt2, nomatch=0], times=100, unit="ms") # fusion de data tables et comparaison avec dplyr (code absent du livre) nrOfRows <- 5*1e7 set.seed(123) x1 <- data.table( num = 1:nrOfRows, flag = sample(c(TRUE, FALSE, NA), nrOfRows, replace = TRUE), etat1 = factor(sample(state.name, nrOfRows, replace = TRUE)), mt1 = runif(nrOfRows, 0.0, 100), stringsAsFactors = FALSE, key = "num") head(x1,10) set.seed(456) x2 <- data.table( num = 1:nrOfRows, etat2 = factor(sample(state.name, nrOfRows, replace = TRUE)), mt2 = runif(nrOfRows, 0.0, 100), stringsAsFactors = FALSE, key = "num") head(x2,10) system.time(x <- merge(x1, x2)) # observations communes head(x) rm(x1);rm(x2);rm(x) library(dplyr) set.seed(123) y1 <- data.frame( num = 1:nrOfRows, flag = sample(c(TRUE, FALSE, NA), nrOfRows, replace = TRUE), etat1 = factor(sample(state.name, nrOfRows, replace = TRUE)), mt1 = runif(nrOfRows, 0.0, 100), stringsAsFactors = FALSE) head(y1,10) set.seed(456) y2 <- data.frame( num = 1:nrOfRows, etat2 = factor(sample(state.name, nrOfRows, replace = TRUE)), mt2 = runif(nrOfRows, 0.0, 100), stringsAsFactors = FALSE) head(y2,10) system.time(y <- y1 %>% left_join(y2, by = c("num"))) head(y) rm(y1);rm(y2);rm(y) # "si alors" rapide de data.table et comparaison avec dplyr (code absent du livre) nrOfRows <- 1e6 set.seed(123) x <- data.table( num = 1:nrOfRows, orig = factor(sample(state.name, nrOfRows, replace = TRUE)), dest = factor(sample(state.name, nrOfRows, replace = TRUE)), stringsAsFactors = FALSE) system.time(x[, c("Wyoming_to_Tennessee","Indiana_to_Iowa") := list(fifelse(x$orig=="Wyoming" & x$dest=="Tennessee",1,0),fifelse(x$orig=="Indiana" & x$dest=="Iowa",1,0))]) system.time(x[, c("Wyoming_to_Tennessee","Indiana_to_Iowa") := list(ifelse(x$orig=="Wyoming" & x$dest=="Tennessee",1,0),ifelse(x$orig=="Indiana" & x$dest=="Iowa",1,0)),]) system.time(x[, c("Wyoming_to_Tennessee","Indiana_to_Iowa") := list((x$orig=="Wyoming" & x$dest=="Tennessee")*1L,(x$orig=="Indiana" & x$dest=="Iowa")*1L)]) head(x) set.seed(123) x <- data.frame( num = 1:nrOfRows, orig = factor(sample(state.name, nrOfRows, replace = TRUE)), dest = factor(sample(state.name, nrOfRows, replace = TRUE)), stringsAsFactors = FALSE) system.time(x <- x %>% mutate(Wyoming_to_Tennessee = if_else(x$orig=="Wyoming" & x$dest=="Tennessee",1,0), Indiana_to_Iowa = if_else(x$orig=="Indiana" & x$dest=="Iowa",1,0))) head(x) # calcul d'une nouvelle variable : nombre de vols de l'avion dans les dernières 24 heures system.time(v2008.dt <- fread("D:/Data/vols/2008.csv", nrows=1000000)) v2008.dt <- v2008.dt[TailNum > 0] nrow(v2008.dt) # format date library(fasttime) library(stringr) system.time(v2008.dt[,':='(Temps = fastPOSIXct(paste(paste0(Year,"-",str_pad(Month,2,pad = "0"),"-",str_pad(DayofMonth,2,pad = "0")),paste0(substr(str_pad(CRSDepTime,4,pad = "0"),1,2),":",substr(str_pad(CRSDepTime,4,pad = "0"),3,4),":","00")),"GMT"))]) system.time(setkey(v2008.dt, TailNum, Temps)) laps <- 1 # 1 jour : laps de temps pour le comptage du nombre de vols nb_vols_max <- 30 # calcul rapide de nb_vols system.time(v2008.dt[, paste0("num_lag", 1:nb_vols_max) := lapply(1:nb_vols_max, function(i) shift(TailNum, i, type = "lag"))] [, paste0("time_lag", 1:nb_vols_max) := lapply(1:nb_vols_max, function(i) shift(Temps, i, type = "lag"))] [, paste0("time_diff", 1:nb_vols_max) := lapply(1:nb_vols_max, function(i) difftime(Temps, get(paste0("time_lag", i)), units = "days"))] [, nb_vols := 1 + eval(parse(text=paste(paste0("ifelse(is.na(time_diff",1:nb_vols_max,")|(num_lag",1:nb_vols_max,"!= TailNum),0,(time_diff",1:nb_vols_max," <= laps))"),collapse="+")))] [, paste0("num_lag", 1:nb_vols_max) := NULL] [, paste0("time_diff", 1:nb_vols_max) := NULL] [, paste0("time_lag", 1:nb_vols_max) := NULL]) table(v2008.dt$nb_vols) # calcul de nb_vols avec une boucle beaucoup plus longue nb_vols <- rep(1, nrow(v2008.dt)) #system.time(setkey(v2008.dt, TailNum, Temps)) system.time( for (i in 2:nrow(v2008.dt)) { j <- 1 while ((j < i) && (difftime(v2008.dt$Temps[i], v2008.dt$Temps[i-j], units="days") <= laps) && (v2008.dt$TailNum[i]==v2008.dt$TailNum[i-j])) { nb_vols[i] <- nb_vols[i] + 1 j <- j + 1 } } ) table(nb_vols) identical(nb_vols, v2008.dt$nb_vols) all.equal(nb_vols, v2008.dt$nb_vols, check.attributes = F) # transformation de cas en variables (passage en format large) (code absent du livre) (dt1 <- data.table(A = letters[1:5], no = 1:2, dateA = sample(2011:2020), dateB = sample(2011:2020), key = "A,no")) (dt2 <- dcast(dt1, A ~ no, value.var = c("dateA","dateB"))) # transformation de variables en cas (passage en format long) (code absent du livre) (dt3 <- melt(dt2, id.vars = c("A"), measure.vars = patterns("^dateA", "^dateB"), variable.name = "no", value.name = c("dateA", "dateB"))) # Section 8.8.4 : Lecture et écriture rapide avec le package fst # --------------------------------------------------------------------------------------------------------- # jeu de données avec 1 million de lignes avec des variables de types différents nrOfRows <- 1e6 x <- data.frame( Integers = 1:nrOfRows, # integer Logicals = sample(c(TRUE, FALSE, NA), nrOfRows, replace = TRUE), # logical Text = factor(sample(state.name, nrOfRows, replace = TRUE)), # text Numericals = runif(nrOfRows, 0.0, 100), # numericals stringsAsFactors = FALSE) head(x) summary(x) #install.packages('fst') library(fst) library(data.table) # écriture system.time(write.fst(x, "dataset.fst")) system.time(fwrite(x, "dataset.dt")) system.time(write.csv(x, "dataset.csv", row.names = FALSE)) # lecture system.time(y <- read.fst("dataset.fst")) dim(y) identical(x,y) all.equal(x, y, check.attributes = FALSE) system.time(yt <- fread("dataset.dt", stringsAsFactors = TRUE)) dim(yt) all.equal(x, yt, check.attributes = FALSE) system.time(yc <- read.csv("dataset.csv", stringsAsFactors = TRUE)) dim(yc) all.equal(x, yc, check.attributes = FALSE) # compression maximale write.fst(x, "dataset.fst.max", 100) # compression minimale write.fst(x, "dataset.fst.min", 0) # Section 8.8.6 : Traitements sans charger toutes les données en mémoire # --------------------------------------------------------------------------------------------------------- # lecture avec package sqldf #install.packages('sqldf') library(sqldf) f <- file("D:/Data/vols/2008.csv") v2008.sql <- read.csv.sql(sql = "select * from f", eol="\n") v2008.sql <- read.csv.sql(sql = "select * from f where UniqueCarrier = 'AA' and Origin='JFK'", eol="\n") class(v2008.sql) dim(v2008.sql) # lecture avec package LaF #install.packages('LaF') library(LaF) datamodel <- detect_dm_csv("D:/Data/vols/2008.csv", sep=",", header=TRUE) datamodel$columns$type # type des colonnes du fichier csv v2008.laf <- laf_open_csv("D:/Data/vols/2008.csv", column_types=datamodel$columns$type, skip=1) v2008.laf <- laf_open(datamodel) # on peut utiliser directement laf_open si on a créé un "model" par detect_dm_csv v2008.lf <- v2008.laf[1:100,1:23] # on ne peut extraire que les 23 premières colonnes car les suivantes contiennent des valeurs manquantes v2008.lf <- v2008.laf[,] # Section 8.8.7 : Traitements sur le disque dur – le package ff # --------------------------------------------------------------------------------------------------------- #install.packages('ff') #install.packages('ffbase') library(ff) getOption("fftempdir") # répertoire temporaire options(fftempdir = "D:/Data/ff") # choix du répertoire temporaire x <- integer(1e+08) object.size(x) x <- integer(1e+9) x <- ff(0, length = 1e+9) object.size(x) x <- 1:1000000 object.size(x) mean(x) x_ff <- ff(x) object.size(x_ff) mean(x_ff) library(ffbase) mean(x_ff) # fonction donnant la date du premier vol en nombre de mois birthmonth <- function (y) { yearInds <- which.min(y[,'Year']) minYear <- min(y[,'Year'], na.rm = TRUE) minMonth <- min(y[yearInds,'Month'], na.rm = TRUE ) return (12*minYear + minMonth - 1) } vols <- read.csv.ffdf(file="D:/Data/vols/2008.csv",header=TRUE) vols <- read.csv.ffdf(file="C:/temp/Vols/2008.csv",header=TRUE) library(ffbase) # calcul de l'âge de l'avion au moment du vol time <- proc.time() planeindices <- split(1:nrow(vols), as.ram(vols$TailNum)) planeStart <- sapply(planeindices, function(i) birthmonth(vols[i, c('Year','Month'), drop=FALSE])) planes <- ff(planeStart) vols$Age <- vols$Year*as.integer(12) + vols$Month - planes[vols$TailNum] proc.time() - time # Section 8.9 : Régression sur données massives avec R # --------------------------------------------------------------------------------------------------------- # régressions sur les données : http://stat-computing.org/dataexpo/2009/ system.time(v2008 <- read.csv(file="D:/Data/vols/2008.csv", header=TRUE)) # régression linéaire sur un data frame de 7 millions de lignes system.time(lin <- lm(ArrDelay ~ Distance, data=v2008)) #system.time(lin <- glm(ArrDelay ~ Distance, data=v2008)) summary(lin) # régression avec biglm sur data frame library(biglm) system.time(biglin <- biglm(ArrDelay ~ Distance, data=v2008)) #system.time(biglin <- bigglm(ArrDelay ~ Distance, data=v2008, family=gaussian(), chunksize=1000000)) summary(biglin) # création d'un fichier ff library(ff) v2008.ff <- read.csv.ffdf(file="D:/Data/vols/2008.csv",header=TRUE) object.size(v2008.ff) object.size(v2008) # régression linéaire sur un data frame ffdf library(biglm) system.time(fflin <- biglm(ArrDelay ~ Distance, data=v2008.ff)) summary(fflin) # dépassement de la capacité mémoire system.time(v2007 <- read.csv(file="D:/Data/vols/2007.csv", header=TRUE)) system.time(v20072008 <- rbind(v2007, v2008)) # traitement avec package ff v2007.ff <- read.csv.ffdf(file="D:/Data/vols/2007.csv",header=TRUE) library(ffbase) system.time(vols <- ffdfappend(v2008.ff, v2007.ff)) class(vols) nrow(vols) avions <- as.data.frame(vols) system.time(fflin <- bigglm(ArrDelay ~ Distance, data=vols, chunksize=100000)) system.time(fflin <- biglm(ArrDelay ~ Distance, data=vols)) summary(fflin) # régression avec biglm et update library(data.table) system.time(v2008.dt <- fread("D:/Data/vols/2008.csv")) system.time(biglin <- bigglm(ArrDelay ~ Distance, data=v2008.dt, family=gaussian(), chunksize=1000000)) system.time(biglin <- biglm(ArrDelay ~ Distance, data=v2008.dt)) summary(biglin) rm(v2008.dt) gc() system.time(v2007.dt <- fread("D:/Data/vols/2007.csv")) system.time(biglin <- update(biglin, v2007.dt)) summary(biglin) # régression linéaire sur 5 prédicteurs system.time(fflin <- bigglm(ArrDelay ~ DepDelay + Distance + Year + ActualElapsedTime + UniqueCarrier, data=vols, chunksize=100000)) system.time(fflin <- biglm(ArrDelay ~ DepDelay + Distance + Year + ActualElapsedTime + UniqueCarrier, data=vols)) summary(fflin) # concaténation de plusieurs années v2006.ff <- read.csv.ffdf(file="D:/Data/vols/2006.csv",header=TRUE) v2005.ff <- read.csv.ffdf(file="D:/Data/vols/2005.csv",header=TRUE) v2004.ff <- read.csv.ffdf(file="D:/Data/vols/2004.csv",header=TRUE) time <- proc.time() vols <- ffdfappend(v2008.ff, v2007.ff) vols <- ffdfappend(vols,v2006.ff) vols <- ffdfappend(vols,v2005.ff) vols <- ffdfappend(vols,v2004.ff) proc.time() - time nrow(vols) system.time(fflin <- bigglm(ArrDelay ~ DepDelay + Distance + Year + ActualElapsedTime + UniqueCarrier, data=vols, chunksize=100000)) summary(fflin) system.time(fflin <- biglm(ArrDelay ~ DepDelay + Distance + Year + ActualElapsedTime + UniqueCarrier, data=vols)) summary(fflin) # avec H2O library(h2o) localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, max_mem_size = '4g', ice_root="D:/Data/h2o", nthreads=2) system.time(vols <- h2o.importFolder(path = "D:/Data/volszip", header=TRUE)) nrow(vols) system.time(lmh2o <- h2o.glm(y = "ArrDelay", x = c("DepDelay", "Distance", "Year", "ActualElapsedTime","UniqueCarrier"), training_frame = vols, family = "gaussian", nfolds = 0, alpha = 0, lambda = 0)) print(lmh2o) # Section 8.10.1 : Package biglm sur data frame # --------------------------------------------------------------------------------------------------------- memory.size(max=NA) memory.limit(size=1024) set.seed(123) matest <- matrix(rnorm(30000000,0,2),10000000,3) colnames(matest) <- c("v1","v2","v3") object.size(matest) # 240 Mo # modèle linéaire avec lm system.time(lin <- lm(v1 ~ v2 + v3, data=as.data.frame(matest))) summary(lin) # modèle linéaire avec bigglm "de base" (en RAM) library(biglm) # modèle linéaire : méthode des MCO (calcul matriciel) system.time(biglin <- biglm(v1 ~ v2 + v3, data=as.data.frame(matest))) summary(biglin) # modèle linéaire généralisé : méthode de Fisher de l’algorithme de Newton-Raphson system.time(biglin <- bigglm(v1 ~ v2 + v3, data=as.data.frame(matest), family=gaussian(), chunksize=10000)) system.time(biglin <- bigglm(v1 ~ v2 + v3, data=as.data.frame(matest), family=gaussian(), chunksize=100000)) summary(biglin) # Section 8.10.2 : Package biglm sur big matrix # --------------------------------------------------------------------------------------------------------- # modèle linéaire avec bigglm sur big matrix install.packages('bigmemory.sri') install.packages("bigmemory") install.packages("biganalytics") install.packages("bigalgebra") library(bigmemory) library(biganalytics) bigtest <- as.big.matrix(matest) system.time(bigmat <- biglm.big.matrix(v1 ~ v2 + v3, data=bigtest)) # KO sur Windows (du moins jusqu'en juin 2020) : # KO : Error in eval(expr, envir, enclos) : objet 'v1' introuvable system.time(bigmat <- bigglm.big.matrix(v1 ~ v2 + v3, data=bigtest, chunksize=100000)) # OK avec chunksize summary(bigmat) system.time(bigmat <- bigglm.big.matrix(v1 ~ v2 + v3, data=bigtest, family=gaussian(), chunksize=100000)) # OK avec chunksize summary(bigmat) # Section 8.10.3 : Packages ff et ffbase sur data frame ffdf # --------------------------------------------------------------------------------------------------------- # modèle linéaire avec packages ff, ffbase et biglm library(ff) options(fftempdir = "D:/Data/ff") # choix du répertoire temporaire ffdftest <- as.ffdf(as.data.frame(matest)) # création d'un objet ffdf à partir d'un data frame library(biglm) system.time(fflin <- biglm(v1 ~ v2 + v3, data=ffdftest)) summary(fflin) library(ffbase) system.time(fflin <- bigglm(v1 ~ v2 + v3, data=ffdftest)) system.time(fflin <- bigglm(v1 ~ v2 + v3, data=ffdftest, chunksize=100000)) summary(fflin) # modèle linéaire avec packages biglm et ff (données trop grosses pour la mémoire vive) set.seed(123) matest <- matrix(rnorm(300000000,0,2),100000000,3) fftest <- ff(rnorm(300000000,0,2),dim=c(100000000,3)) # régression biglm sur objet ff : impossible sur un objet ff matrix => il faut créer un objet ffdf ffdftest <- as.ffdf(fftest) colnames(ffdftest) <- c("v1","v2","v3") head(ffdftest) system.time(fflin <- bigglm(v1 ~ v2 + v3, data=ffdftest, chunksize=100000)) summary(fflin) system.time(fflin <- biglm(v1 ~ v2 + v3, data=ffdftest)) # mécanisme chunkwise pour biglm formule <- v1 ~ v2 + v3 time <- proc.time() for (i in chunk(ffdftest, by=100000)){ if (i[1]==1) { biglmfit <- biglm(formule, data=ffdftest[i,]) } else { biglmfit <- update(biglmfit, ffdftest[i,]) } } proc.time() - time summary(biglmfit) # Section 8.10.4 : Package speedglm sur data frame # --------------------------------------------------------------------------------------------------------- matest <- matrix(rnorm(30000000,0,2),10000000,3) # 240 Mo colnames(matest) <- c("v1","v2","v3") library(speedglm) system.time(speedlin <- speedlm(v1 ~ v2 + v3, data=as.data.frame(matest))) system.time(speedlin <- speedglm(v1 ~ v2 + v3, data=as.data.frame(matest),family=gaussian())) system.time(speedlin <- speedglm(v1 ~ v2 + v3, data=as.data.frame(matest))) summary(speedlin) # Section 8.10.6 : Comparaison des temps de calcul d'une régression linéaire # --------------------------------------------------------------------------------------------------------- library(bigmemory) library(biganalytics) library(speedglm) library(RcppArmadillo) library(microbenchmark) set.seed(123) matest <- matrix(rnorm(3000000,0,2),1000000,3) # 24 Mo matest <- matrix(rnorm(15000000,0,2),5000000,3) # 120 Mo matest <- matrix(rnorm(30000000,0,2),10000000,3) # 240 Mo colnames(matest) <- c("v1","v2","v3") test <- as.data.frame(matest) bigtest <- as.big.matrix(matest) compare <- microbenchmark(lm(v1 ~ v2 + v3, data=test), glm(v1 ~ v2 + v3, data=test), biglm(v1 ~ v2 + v3, data=test), bigglm(v1 ~ v2 + v3, data=test, family=gaussian()), speedlm(v1 ~ v2 + v3, data=test), speedglm(v1 ~ v2 + v3, data=test, family=gaussian()), fastLm(v1 ~ v2 + v3, data=test, family=gaussian()), biglm.big.matrix(v1 ~ v2 + v3, data=bigtest, chunksize=100000), bigglm.big.matrix(v1 ~ v2 + v3, data=bigtest, chunksize=100000), times=30, unit="ms") compare # Section 8.12 : Traitement des matrices creuses # --------------------------------------------------------------------------------------------------------- m1 <- matrix(0, nrow = 10000, ncol = 10000) object.size(m1) /(1024 ^ 2) library('Matrix') m2 <- Matrix(0, nrow = 10000, ncol = 10000, sparse = TRUE) object.size(m2) /(1024 ^ 2) set.seed(235) system.time(m2[round(10000*runif(1000,0,1)),round(10000*runif(1000,0,1))] <- 1) system.time(m2[sample(1e4)[1:1e3], sample(1e4)[1:1e3]] <- 1) nnzero(m2) / (dim(m2)[1]*dim(m2)[2]) object.size(m2) /(1024 ^ 2) m2[1:30,1:30] image(m2[1:30,1:30]) # création d'une matrice sparse en indiquant les coordonnées i,j des éléments non nuls # les valeurs x de ces éléments non nuls, et la dimension de la matrice i <- sample(1e4,1e6,replace=T) j <- sample(1e4,1e6,replace=T) x <- 1 x <- rnorm(1e6) x <- rpois(1e6,100) system.time(m3 <- sparseMatrix(i=i, j=j, x=x, dims=list(10000,10000))) nnzero(m3) / (dim(m3)[1]*dim(m3)[2]) object.size(m3) /(1024 ^ 2) m3[1:30,1:30] image(m3[1:30,1:30]) # création des indicatrices correspondant aux niveaux d’un facteur credit = read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data", h=FALSE) credit$V21 <- as.factor(credit$V21) library(glmnet) y <- credit[,"V21"] x <- model.matrix( ~ . -1 , data=credit[,-which(names(credit)=="V21")]) object.size(x) x <- sparse.model.matrix( ~ . -1 , data=credit[,-which(names(credit)=="V21")]) object.size(x) set.seed(235) cvfit <- cv.glmnet(x, y, alpha=0.001, family = "binomial", type="auc", nlambda=100) cvfit # --------------------------------------------------------------------------------------------------------- # Chapitre 9 : Big Data avec d’autres logiciels # --------------------------------------------------------------------------------------------------------- # Section 9.2 : Utilisation d’une bibliothèque optimisée de calcul # --------------------------------------------------------------------------------------------------------- # comparaison entre R CRAN et R Open A <- matrix(rnorm(10000^2),10000) dim(A) system.time(A%*%A) system.time(solve(A)) system.time(qr.coef(qr(A, LAPACK = FALSE),diag(1,nrow(A)))) # paramètre par défaut system.time(qr.coef(qr(A, LAPACK = TRUE),diag(1,nrow(A)))) # paramètre plus rapide # comparaison (plus complète) des méthodes pour la régression library(bigmemory) library(biganalytics) library(speedglm) library(RcppArmadillo) library(ff) library(ffbase) library(microbenchmark) set.seed(123) m <- 100000 n <- 100 matest <- matrix(rnorm(m*n,0,2), m, n) object.size(matest)/1024^2 # 76 Mo dim(matest) y <- rnorm(m) # variable à expliquer x <- cbind(1, matest) # vecteur de variables explicatives (y compris la constante) test <- data.frame(matest, y = y) # lm nécessite un data frame et non une matrice # remarque : biglm n'admet pas de formule de la forme y ~., et il faut donc écrire la liste des prédicteurs nm <- names(test) f <- as.formula(paste("y ~", paste(nm[!nm %in% "y"], collapse = " + "))) bigtest <- as.big.matrix(cbind(matest, y)) options(bigmemory.allow.dimnames=TRUE) colnames(bigtest) <- names(test) # package ff options(fftempdir = "D:/Data/ff") # choix du répertoire temporaire ffdftest <- as.ffdf(test) # nombre de threads utilisés par MKL sur R Open getMKLthreads() setMKLthreads(2) # comparaison des temps de calcul compare <- microbenchmark(lm(y ~ ., data=test), lm.fit(x,y), biglm(f, data=test), bigglm(f, data=test, family=gaussian()), speedlm(y ~ ., data=test), speedglm(y ~ ., data=test, family=gaussian()), fastLm(y ~ ., data=test, family=gaussian()), biglm.big.matrix(f, data=bigtest, chunksize=100000), bigglm.big.matrix(f, data=bigtest, chunksize=100000), biglm(f, data=ffdftest), bigglm(f, data=ffdftest, chunksize=100000), solve(crossprod(x), crossprod(x,y)), solve(t(x)%*%x, t(x)%*%y), qr.coef(qr(x, LAPACK = TRUE), y), qr.coef(qr(x, LAPACK = FALSE), y), times=100, unit="ms") compare # Section 9.4.2 : Installation du package h2o # --------------------------------------------------------------------------------------------------------- # installation à partir du CRAN install.packages('h2o') # installation de la version la plus récente # suppression anciennes versions if ("package:h2o" %in% search()) { detach("package:h2o", unload=TRUE) } if ("h2o" %in% rownames(installed.packages())) { remove.packages("h2o") } # dépendances pkgs <- c("RCurl","jsonlite") for (pkg in pkgs) { if (! (pkg %in% rownames(installed.packages()))) { install.packages(pkg) } } # version la plus récente install.packages("h2o", type="source", repos=(c("http://h2o-release.s3.amazonaws.com/h2o/latest_stable_R"))) # Sections 9.4.3 et 9.4.4 : Démarrage et lecture de données avec le package h2o # --------------------------------------------------------------------------------------------------------- library(h2o) localH2O <- h2o.init() localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, max_mem_size = '4g', ice_root="D:/Data/h2o", nthreads=2) # lecture d'un fichier system.time(vols <- h2o.importFile(path = "D:/Data/vols/2008.csv")) system.time(vols <- h2o.importFile(path = "D:/Data/vols/2008.csv", header=TRUE, sep=",")) system.time(v2008 <- h2o.importFile(path = "D:/Data/volszip/2008.zip")) system.time(vols <- h2o.importFolder(path = "D:/Data/volszip")) head(vols) classes <- c("Enum", "Numeric", "Numeric", …, "Numeric", "Enum") base <- h2o.importFile(path = "D:\\Data\\base.csv", col.types = classes) h2o.ls() # suppression d'un objet du cluster H2O h2o.rm(vols) # suppression de tous les objets du cluster H2O h2o.rm(as.character(h2o.ls()[,1])) # arrêt du cluster h2o.shutdown() h2o.shutdown(prompt=FALSE) # Section 9.4.5 : Exécution de fonctions dans un cluster h2o # --------------------------------------------------------------------------------------------------------- h2o.table(vols$UniqueCarrier) head(h2o.table(vols$UniqueCarrier), 30) system.time(vols$Late <- (vols$ArrDelay > 15)) head(h2o.table(vols$UniqueCarrier, vols$Late), 10) system.time(cross <- h2o.group_by(vols, by = c("UniqueCarrier","Origin"), nrow("Distance"), mean("Distance"))) head(cross, 10) # Section 9.4.6 : Régression avec h2o # --------------------------------------------------------------------------------------------------------- system.time(vols$Hour <- floor(vols$CRSDepTime/100)) system.time(glmh2o <- h2o.glm(y = "Late", x = c("Distance", "Year", "Hour", "CRSElapsedTime", "UniqueCarrier"), training_frame = vols, family = "binomial", nfolds = 0, alpha = 0, lambda = 0)) print(glmh2o) # Section 9.5.1 : Le package sparklyr # --------------------------------------------------------------------------------------------------------- install.packages('sparklyr') devtools::install_github("rstudio/sparklyr") spark_install(version = "2.4.3") spark_available_versions() library(sparklyr) config <- spark_config() config$spark.executor.cores <- 4 config$spark.executor.memory <- "2G" #config$spark.driver.memory <- "4G" Sys.setenv(SPARK_HOME="D:/Spark/spark-2.4.5-bin-hadoop2.7") Sys.setenv(JAVA_HOME="C:/Program Files/Java/jdk1.8.0_181") config[["spark.r.command"]] <- "C:/Program Files/R/R-3.6.1/bin/Rscript.exe" #sc <- spark_connect(master = "local", config = config, version = "2.2.1") sc <- spark_connect(master = "local", config = config) #sc <- spark_connect(master = "local") # chargement dans Spark à partir d'un fichier CSV system.time(v2008_tbl <- sparklyr::spark_read_csv(sc, name="v2008", path="D:/Data/vols/2008.csv", header = TRUE)) sdf_schema(v2008_tbl) sdf_describe(v2008_tbl) # chargement dans Spark à partir d'un data frame R system.time(v2008 <- read.csv("D:/Data/vols/2008.csv", header=TRUE, stringsAsFactors=TRUE)) system.time(v2008_tbl <- copy_to(sc, v2008, overwrite = TRUE)) # liste des tables src_databases(sc) # Spark SQL DBI::dbGetQuery(sc, "SELECT TailNum, COUNT(TailNum), AVG(Distance) FROM v2008 GROUP BY TailNum LIMIT 10") DBI::dbGetQuery(sc, "SELECT TailNum, COUNT(TailNum), AVG(Distance) FROM v2008 GROUP BY TailNum") library(dplyr) v2008_tbl %>% group_by(TailNum) %>% summarise(n = n(), mean = mean(Distance)) # régression linéaire simple system.time(fit <- ml_linear_regression(x = v2008_tbl, response = "ArrDelay", features = c("Distance"))) summary(fit) # copie d’un data frame dans Spark v2008 <- read.csv("D:/Data/vols/2008.csv", header=TRUE, stringsAsFactors=TRUE) system.time(v2008_tbl <- copy_to(sc, v2008, overwrite = TRUE)) # régression logistique system.time(flights <- v2008_tbl %>% mutate(Late = ifelse(ArrDelay > "15",1,0), Hour = floor(CRSDepTime/100))) system.time(fit <- flights %>% ml_generalized_linear_regression(response = "Late", features = c("Distance", "Hour", "CRSDepTime", "UniqueCarrier"), family = binomial("logit"))) system.time(summary(fit)) # la fonction spark_apply devrait permettre d'exécuter dans un cluster Spark n'importe quelle fonction d'un package de R conçue pour s'appliquer à un data frame # spark_apply() appliquera cette fonction à un objet Spark (Spark DataFrame) même s'il est distribué dans plusieurs partitions d'un cluster # spark_apply() appliquera cette fonction à chaque partition du cluster et renverra le résultat dans un unique Spark DataFrame # on a parfois le message d'erreur : Error in file(con, "r") : impossible d'ouvrir la connexion iris_tbl <- copy_to(sc, iris) spark_apply(iris_tbl, function(e) print(MASS::lda(Species ~ ., e))) spark_apply(iris_tbl, function(e) coefficients(MASS::lda(Species ~ ., e))) spark_apply(iris_tbl, function(e) broom::tidy(lm(Petal_Length ~ Petal_Width, e)), names = c("term", "estimate", "std.error", "statistic", "p.value"), group_by = "Species") spark_apply(iris_tbl, function(e) broom::tidy(glm(Species ~ ., e, family=binomial(link = "logit"))), names = c("term", "estimate", "std.error", "statistic", "p.value")) library(MASS) print(lda(Species ~ ., data=iris)) coefficients(lda(Species ~ ., data=iris)) # déconnexion spark_disconnect(sc) # Section 9.5.2 : Le package SparkR # --------------------------------------------------------------------------------------------------------- # installation SparkR à partir du CRAN install.packages('SparkR') library(SparkR) # cette ligne fait pointer la variable d’environnement SPARK_HOME vers le dossier où Spark est installé Sys.setenv(SPARK_HOME="D:/Spark/spark-2.4.5-bin-hadoop2.7") # cette ligne fait pointer la variable d’environnement JAVA_HOME vers le dossier où Java est installé Sys.setenv(JAVA_HOME="C:/Program Files/Java/jdk1.8.0_181") # création d'un cluster (d'un "SparkContext") sc <- sparkR.session(master = "local", sparkConfig = list(spark.driver.memory = "2g", spark.driver.cores = "4")) # création d'un Spark DataFrame df <- data.frame(name=c("John", "Smith", "Sarah"), age=c(19, 23, 18)) spdf <- createDataFrame(df) printSchema(spdf) # lecture d'un fichier schema <- structType( structField("SepalLength", "double", FALSE), structField("SepalWidth", "double", FALSE), structField("PetalLength", "double", FALSE), structField("PetalWidth", "double", FALSE), structField("Species", "string", FALSE)) write.csv(iris, file = "iris.csv", row.names = F) irisDF <- read.df("iris.csv", "csv", header="true", schema=schema) head(irisDF) # régression linéaire system.time(model <- spark.glm(irisDF, SepalLength ~ SepalWidth + Species, family = "gaussian")) system.time(model <- glm(SepalLength ~ SepalWidth + Species, data = irisDF, family = "gaussian")) summary(model) # traitement de données massives system.time(v2008 <- read.df("D:/Data/vols/2008.csv", "csv", header = "true")) system.time(v2008 <- read.df("D:/Data/vols/2008.csv", "csv", header = "true", inferSchema = "true", na.strings = "NA")) system.time(v2008 <- loadDF("D:/Data/vols/2008.csv", source="csv", header = "true", inferSchema = "true", na.strings = "NA")) system.time(v2008$Late <- ifelse(v2008$ArrDelay > "15",1,0)) printSchema(v2008) head(v2008) # remplacement ou suppression des valeurs manquantes : indispensable pour la variable à expliquer v2008 <- dropna(v2008, cols = "ArrDelay") describe(v2008, "ArrDelay") system.time(model <- spark.glm(v2008, ArrDelay ~ Distance, family = "gaussian")) system.time(model <- spark.glm(v2008, ArrDelay ~ DepDelay + Distance + ActualElapsedTime + UniqueCarrier, family = "gaussian")) time <- proc.time() summary(model) proc.time() - time # arbres de décision system.time(model <- spark.decisionTree(v2008, ArrDelay ~ Distance, type = "regression", maxDepth = 4, maxBins = 16)) system.time(model <- spark.decisionTree(v2008, ArrDelay ~ DepDelay + Distance + ActualElapsedTime + UniqueCarrier, type = "regression", maxDepth = 4, maxBins = 16)) pred.tree <- predict(model, v2008) # gradient boosting system.time(model <- spark.gbt(v2008, ArrDelay ~ Distance, type = "regression", maxDepth = 5, maxBins = 16)) # forêts aléatoires system.time(model <- spark.randomForest(v2008, ArrDelay ~ Distance, type = "regression", maxDepth = 5, maxBins = 16)) # perceptron system.time(model <- spark.mlp(v2008, CRSArrTime ~ Distance, blockSize = 128, layers = c(4,3), solver = "l-bfgs", maxIter = 100, tol = 0.5, stepSize = 1, seed = 123)) # SVM system.time(model <- spark.svmLinear(v2008, Late ~ Distance, regParam = 0.5)) # Interface Web Spark http://localhost:4040/ # arrêt de SparkR sparkR.stop() # Section 9.8 : Comparaison des logiciels # --------------------------------------------------------------------------------------------------------- library(ff) library(ffbase) options(fftempdir = "D:/Data/ff") # lecture optimisée avec package ff entete = read.csv(file="D:/Data/vols/2004.csv", header = TRUE, nrows = 1000) headclasses = sapply(entete, class) # import des fichiers CSV time <- proc.time() v2004.ff <- read.csv.ffdf(file="D:/Data/vols/2004.csv",header=TRUE, colClasses=headclasses, BATCHBYTES = 20*getOption("ffbatchbytes")) v2005.ff <- read.csv.ffdf(file="D:/Data/vols/2005.csv",header=TRUE, colClasses=headclasses, BATCHBYTES = 20*getOption("ffbatchbytes")) v2006.ff <- read.csv.ffdf(file="D:/Data/vols/2006.csv",header=TRUE, colClasses=headclasses, BATCHBYTES = 20*getOption("ffbatchbytes")) v2007.ff <- read.csv.ffdf(file="D:/Data/vols/2007.csv",header=TRUE, colClasses=headclasses, BATCHBYTES = 20*getOption("ffbatchbytes")) v2008.ff <- read.csv.ffdf(file="D:/Data/vols/2008.csv",header=TRUE, colClasses=headclasses, BATCHBYTES = 20*getOption("ffbatchbytes")) proc.time() - time # concaténation de plusieurs années time <- proc.time() vols <- ffdfappend(v2008.ff,v2007.ff) vols <- ffdfappend(vols,v2006.ff) vols <- ffdfappend(vols,v2005.ff) vols <- ffdfappend(vols,v2004.ff) proc.time() - time nrow(vols) # régression linéaire library(biglm) system.time(fflin <- bigglm(ArrDelay ~ DepDelay+Distance+Year+ActualElapsedTime+UniqueCarrier, data=vols, chunksize=1000000)) summary(fflin) # régression logistique time <- proc.time() vols$Hour <- floor(vols$CRSDepTime/100) vols$Late <- (vols$ArrDelay > 15) proc.time() - time system.time(fflin <- bigglm(Late ~ Distance+Year+Hour+CRSElapsedTime+UniqueCarrier, data=vols, family = binomial(), chunksize=1000000, tolerance=1e-4)) summary(fflin) library(h2o) localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, max_mem_size = '4g', ice_root="D:/Data/h2o", nthreads=-1) system.time(vols <- h2o.importFolder(path = "D:/Data/volszip")) system.time(glmh2o <- h2o.glm(y = "ArrDelay", x = c("DepDelay", "Distance", "Year", "ActualElapsedTime", "UniqueCarrier"), training_frame = vols, family = "gaussian", nfolds = 0, alpha = 0, lambda = 0)) print(glmh2o) time <- proc.time() vols$Hour <- floor(vols$CRSDepTime/100) vols$Late <- (vols$ArrDelay > 15) proc.time() - time system.time(glmh2o <- h2o.glm(y = "Late", x = c("Distance", "Year", "Hour", "CRSElapsedTime", "UniqueCarrier"), training_frame = vols, family = "binomial", nfolds = 0, alpha = 0, lambda = 0)) print(glmh2o) # --------------------------------------------------------------------------------------------------------- # Chapitre 10 : Traitement du langage naturel # --------------------------------------------------------------------------------------------------------- # Section 10.3.2 : Identification de la langue # --------------------------------------------------------------------------------------------------------- library(textcat) textcat(c("This is an English sentence.", "Das ist ein deutscher Satz.", "Esta es una frase en español.", "C'est une phrase en français.", "E' una frase in italiano.")) library(franc) sapply(c("This is an English sentence.", "Das ist ein deutscher Satz.", "Esta es una frase en español.", "C'est une phrase en français.", "E' una frase in italiano."), franc) # loi de Zipf du nombre de locuteurs locuteurs <- speakers sapply(1:20, function(n) (locuteurs$speakers[1]/locuteurs$speakers[n+1])) # Section 10.3.3 : Segmentation (tokenisation) # --------------------------------------------------------------------------------------------------------- library(quanteda) tokens("My taylor is rich", what="word") tokens("M. Taylor is rich", what="sentence") texte <- "To be or not to be, that is the question." tokens(texte) tokens_ngrams(tokens(texte), n = 2:3) tokens(texte, remove_punct = TRUE, ngrams = 2, concatenator = " ") tokens(texte, remove_punct = TRUE, ngrams = 2, concatenator = " ", skip = 0:1) tokens("My taylor is rich", remove_punct = TRUE, ngrams = 1:2) # Section 10.3.4 : Identification des catégories grammaticales (étiquetage grammatical) # --------------------------------------------------------------------------------------------------------- # Voir section 10.3.7 # Sections 10.3.5 et 10.3.6 : Entités nommées # --------------------------------------------------------------------------------------------------------- #install.packages('cleanNLP') library(cleanNLP) # Incipit Pride and Prejudice txt <- c(" It is a truth universally acknowledged, that a single man in possession of a good fortune, must be in want of a wife. However little known the feelings or views of such a man may be on his first entering a neighbourhood, this truth is so well fixed in the minds of the surrounding families, that he is considered the rightful property of some one or other of their daughters. 'My dear Mr. Bennet,' said his lady to him one day, 'have you heard that Netherfield Park is let at last?' Mr. Bennet replied that he had not. 'But it is,' returned she; 'for Mrs. Long has just been here, and she told me all about it.' Mr. Bennet made no answer. 'Do you not want to know who has taken it?' cried his wife impatiently. 'You want to tell me, and I have no objection to hearing it.'") # Incipit Moby Dick txt <- c("Call me Ishmael. Some years ago—never mind how long precisely—having little or no money in my purse, and nothing particular to interest me on shore, I thought I would sail about a little and see the watery part of the world. It is a way I have of driving off the spleen and regulating the circulation. Whenever I find myself growing grim about the mouth; whenever it is a damp, drizzly November in my soul; whenever I find myself involuntarily pausing before coffin warehouses, and bringing up the rear of every funeral I meet; and especially whenever my hypos get such an upper hand of me, that it requires a strong moral principle to prevent me from deliberately stepping into the street, and methodically knocking people’s hats off—then, I account it high time to get to sea as soon as I can. This is my substitute for pistol and ball. With a philosophical flourish Cato throws himself upon his sword; I quietly take to the ship. There is nothing surprising in this. If they but knew it, almost all men in their degree, some time or other, cherish very nearly the same feelings towards the ocean with me.", "There now is your insular city of the Manhattoes, belted round by wharves as Indian isles by coral reefs—commerce surrounds it with her surf. Right and left, the streets take you waterward. Its extreme downtown is the battery, where that noble mole is washed by waves, and cooled by breezes, which a few hours previous were out of sight of land. Look at the crowds of water-gazers there.", "Circumambulate the city of a dreamy Sabbath afternoon. Go from Corlears Hook to Coenties Slip, and from thence, by Whitehall, northward. What do you see?—Posted like silent sentinels all around the town, stand thousands upon thousands of mortal men fixed in ocean reveries. Some leaning against the spiles; some seated upon the pier-heads; some looking over the bulwarks of ships from China; some high aloft in the rigging, as if striving to get a still better seaward peep. But these are all landsmen; of week days pent up in lath and plaster—tied to counters, nailed to benches, clinched to desks. How then is this? Are the green fields gone? What do they here?") # Incipit Le Rouge et le Noir txt <- c("La petite ville de Verrières peut passer pour l'une des plus jolies de la Franche-Comté. Ses maisons blanches avec leurs toits pointus de tuiles rouges s'étendent sur la pente d'une colline, dont des touffes de vigoureux châtaigniers marquent les moindres sinuosités. Le Doubs coule à quelques centaines de pieds au-dessous de ses fortifications, bâties jadis par les Espagnols, et maintenant ruinées. Verrières est abritée du côté du nord par une haute montagne, c'est une des branches du Jura. Les cimes brisées du Verra se couvrent de neige dès les premiers froids d'octobre. Un torrent, qui se précipite de la montagne, traverse Verrières avant de se jeter dans le Doubs, et donne le mouvement à un grand nombre de scies à bois; c'est une industrie fort simple et qui procure un certain bien-être à la majeure partie des habitants plus paysans que bourgeois. Ce ne sont pas cependant les scies à bois qui ont enrichi cette petite ville. C'est à la fabrique des toiles peintes, dites de Mulhouse, que l'on doit l'aisance générale qui, depuis la chute de Napoléon, a fait rebâtir les façades de presque toutes les maisons de Verrières.") # Autres exemples txt <- c("Bob Smith gave Alice his login information") txt <- c("Longtemps, je me suis couché de bonne heure.") # lancer l'une des initialisations suivantes avant toute annotation cnlp_init_tokenizers() cnlp_init_spacy() cnlp_init_udpipe() cnlp_init_corenlp() # backend spacy => version > 3.4 de R : installation des fichiers spaCy par la commande suivante # à noter que l'installation du package spacyr ne suffit pas à pouvoir utiliser le backend spaCy avec cleanNLP cnlp_download_spacy(model_name = "en") # backend udpipe => nécessite l'installation du package udpipe et n'extrait ni les entités nommées ni les coréférences (mais les dépendances) install.packages("udpipe") library(udpipe) cnlp_init_udpipe() # backend coreNLP => si pb de téléchargement avec instruction suivante, télécharger directement sur le site : https://stanfordnlp.github.io/CoreNLP/download.html # peut-être préférable à l'instruction cnlp_download_corenlp() qui nécessite d'installer le package coreNLP et de gros fichiers cnlp_download_corenlp() cnlp_download_corenlp(lang="en") cnlp_init_corenlp("en", anno_level = 3, lib_location = NULL, mem = "4g", verbose = FALSE) # analyses en français avec coreNLP # il faut avant le cnlp_init_corenlp() télécharger à l'adresse https://stanfordnlp.github.io/CoreNLP/history.html le fichier Java # stanford-french-corenlp-2016-10-31-models.jar correspondant au français # on trouve aussi des fichiers pour les langues : arabic, chinese , english , german , spanish cnlp_init_corenlp("fr", anno_level = 3, lib_location = NULL, mem = "4g", verbose = FALSE) # annotation annotation <- cnlp_annotate(txt) cnlp_get_token(annotation) # mots du texte avec nature grammaticale et forme lemmatisée cnlp_get_entity(annotation) # entités nommées (annot <- merge(cnlp_get_token(annotation), cnlp_get_entity(annotation), by.x = c("id","tid"), by.y = c("id","tid"), all.x = TRUE)) annot[,c("word","lemma","upos","pos","entity_type","entity")] cnlp_get_dependency(annotation) # dépendances cnlp_get_sentence(annotation) # sentiment de 0 (le plus négatif) à 4 (le plus positif) cnlp_get_coreference(annotation) # coréférences cnlp_get_document(annotation) # métadonnées du document (vector <- cnlp_get_vector(annotation)) # spacy : plongement de mots en 300 dimensions # PCA et topic modeling # Get principal components from the non-proper noun lemmas library(dplyr) tfidf <- cnlp_get_token(annotation) %>% filter(upos %in% c("NOUN", "ADJ", "VERB")) %>% cnlp_get_tfidf() pca_doc <- cnlp_pca(tfidf) # Plot texts using the first two principal components plot(pca_doc$PC1, pca_doc$PC2, col = "white") text(pca_doc$PC1, pca_doc$PC2) # package spacyr # semble préférable à coreNLP pour la lemmatisation en français # https://github.com/quanteda/spacyr #devtools::install_github("quanteda/spacyr", build_vignettes = FALSE) install.packages('spacyr') library(spacyr) spacy_install() # installe spaCy : à lancer avant la première utilisation avec les droits d'administrateur spacy_download_langmodel("fr") # à lancer avant la première utilisation d'une langue autre que l'anglais, avec les droits d'administrateur spacy_initialize() # à lancer avant chaque utilisation spacy_initialize(model = "fr") # à lancer avant chaque utilisation d'une langue autre que l'anglais txt <- "Mr. Smith of moved to San Francisco in December." txt <- "Longtemps, je me suis couché de bonne heure." parsed <- spacy_parse(txt, entity = TRUE) parsed # POS entity_extract(parsed) entity_extract(parsed, type = "all") # entités nommées entity_consolidate(parsed) spacy_finalize() # pour stopper le processus Python # Section 10.3.7 : Lemmatisation # --------------------------------------------------------------------------------------------------------- # https://reaktanz.de/R/pckg/koRpus/koRpus_vignette.html # Le package koRpus peut mener des analyses en anglais, allemand, espagnol, français, italien, néerlandais, portugais et russe, # chacune de ces langues ayant son propre package associé, appelé koRpus.lang.xx, où xx indique la langue # liste des packages disponibles available.koRpus.lang() # installation du package allemand install.koRpus.lang("de") # autre installation possible install.packages("koRpus.lang.de", repos="https://undocumeantit.github.io/repos/l10n/") library(koRpus) available.koRpus.lang() install.koRpus.lang(c("en", "de", "es", "fr", "it")) library(koRpus.lang.fr) tagged.text <- treetag("Longtemps, je me suis couché de bonne heure.", format="obj", treetagger="manual", lang="fr", add.desc = TRUE, TT.options=list(path="C:/TreeTagger/", preset="fr")) tagged.text <- treetag(proust[[1]], format="obj", treetagger="manual", lang="fr", TT.options=list(path="C:/TreeTagger/", preset="fr")) # Un objet de texte tagué contient 5 "slots" (propriétés en langage orienté objet): lang, TT.res (qui contient le text tagué) et desc. # TT.res est un data frame qui contient 8 colonnes (token, tag, lemma, lttr, etc.) et N lignes (une pour chaque token). # Dans R, pour accéder aux slots d'une classe de type S4, on peut utiliser la fonction slot(slot,obj) ou plus simplement utiliser la notation @ txt <- "nous portions des portions d'avocats à des avocats" txt <- "nous sommes au cinéma" txt <- "je suis au cinéma" txt <- "ces cuisiniers excellent dans leur métier" txt <- "la tumeur est détectable par radiographie" tagged.text <- treetag(txt, format="obj", treetagger="manual", lang="fr", TT.options=list(path="C:/TreeTagger/", preset="fr")) tagged.text <- treetag(txt, format="obj", treetagger="manual", lang="fr", TT.options=list(path="C:/TreeTagger/", preset="fr", no.unknown=TRUE)) tagged.text[,-1] tagged.text@desc describe(tagged.text) tagged.text@TT.res taggedText(tagged.text) paste(taggedText(tagged.text)$lemma, collapse=" ") # alternative (très limitée !) à la fonction "treetag" tagged.text <- tokenize("test_lemmatisation.txt", lang="fr") taggedText(tagged.text) # # description du texte summary(tagged.text) plot(tagged.text) # lisibilité ARI(tagged.text) coleman.liau(tagged.text) SMOG(tagged.text) FOG(tagged.text) flesch.kincaid(tagged.text) flesch(tagged.text) flesch(tagged.text, parameters="fr") system.time(readbl.txt <- readability(tagged.text)) # ensemble de tous les indicateurs de lisibilité summary(readbl.txt) # diversité TTR(tagged.text) R.ld(tagged.text) MTLD(tagged.text) maas(tagged.text) system.time(div.txt <- lex.div(tagged.text)) # ensemble de tous les indicateurs de diversité : très long summary(div.txt) lex.div( tagged.text, measure=c("TTR", "MSTTR", "MATTR","HD-D", "MTLD", "MTLD-MA"), char=c("TTR", "MATTR","HD-D", "MTLD", "MTLD-MA") ) # décomposition en syllabes hyph.txt <- hyphen(tagged.text) hyph.txt readbl.txt <- readability(tagged.text, hyphen=hyph.txt) summary(readbl.txt) # autres langues library(koRpus.lang.en) tagged.text <- treetag(txt, format="obj", treetagger="manual", lang="en", TT.options=list(path="C:/TreeTagger/", preset="en")) # Section 10.3.8 : Racinisation # --------------------------------------------------------------------------------------------------------- library(tm) stemDocument("fraises",language = "french") stemDocument("cherche",language = "french") stemDocument("cherchent",language = "french") # Section 10.3.9 : Simplifications # --------------------------------------------------------------------------------------------------------- install.packages('wordnet') library(wordnet) # il faut installer le programme Windows WordNet-2.1 (ou la version Unix WordNet-3.0) téléchargé sur le site https://wordnet.princeton.edu/download/current-version setDict("C:/Program Files (x86)/WordNet/2.1/dict") synonyms("company", "NOUN") filter <- getTermFilter("ExactMatchFilter", "power", TRUE) terms <- getIndexTerms("NOUN", 5, filter) synsets <- getSynsets(terms[[1]]) sapply(synsets, getWord) # Section 10.4.2 : Analyses sur la matrice « documents x termes » # --------------------------------------------------------------------------------------------------------- # liste d'offres d'emploi craigslist <- read.csv("https://raw.githubusercontent.com/h2oai/sparkling-water/rel-1.6/examples/smalldata/craigslistJobTitles.csv") names(craigslist) <- c("category", "text") craigslist$text <- as.character(craigslist$text) head(craigslist) # création corpus library(quanteda) craigs <- corpus(craigslist) summary(craigs, showmeta = TRUE) # création de la matrice « documents x termes » dtm <- dfm(craigs) topfeatures(dtm, 20) dtm <- dfm(craigs, remove = c(stopwords("en")), remove_punct = TRUE, stem = FALSE, tolower = TRUE) topfeatures(dtm, 20) # extrait de la matrice dtm[1:10, 1:10] # termes les plus liés à un terme fix" sim <- textstat_simil(dtm, dtm[,"consultant"], method = "cosine", margin = "features") lapply(as.list(sim), head, 6) # sparsité dtm presDfm <- dfm_trim(dtm, min_termfreq = 200, min_docfreq = 30) presDfm # AFC # remarque : les fonctions textmodel_xx sont dans les versions < 2 du package quanteda # et se trouvent depuis dans le package quanteda.textmodels install.packages('quanteda.textmodels') library(quanteda.textmodels) AFC <- textmodel_ca(presDfm, nd = 2, sparse = TRUE) library(ca) plot.ca(AFC, what=c("none","all")) # meilleur affichage des étiquettes du graphique library(ggplot2) library(ggrepel) ggplot(as.data.frame(AFC$colcoord), aes(x=AFC$colcoord[,1], y=AFC$colcoord[,2])) + geom_point() + geom_text_repel(aes(label=rownames(AFC$colcoord)), hjust=1, vjust=1, size=3, segment.size = 0) # agrégation des documents selon une catégorie dtm.gp <- dfm(craigs, groups = "category", remove = c(stopwords("en")), remove_punct = TRUE, stem = FALSE, tolower = TRUE) dtm.gp dtm.gp[, 1:10] textplot_wordcloud(dtm.gp, min_count = 30, random_order = FALSE, comparison = TRUE, max_words=100, color=gray.colors(3)) # classement de textes y <- craigslist$category set.seed(235) s <- sample(1:nrow(craigslist), nrow(craigslist)*0.5) # classement avec classificateur bayésien naïf (classif <- textmodel_nb(dtm[s,], y[s], prior = "termfreq", distribution = "multinomial")) # matrice de confusion et le taux de bon classement sur l’échantillon de test : table(y[-s], predict(classif, dtm[-s,])) sum(y[-s] == predict(classif, dtm[-s,])) / length(y[-s]) # classement avec forêt aléatoire train <- data.frame(dtm[s,], category=craigslist$category[s]) train$document <- NULL # suppression du no du document (inutile) valid <- data.frame(dtm[-s,], category=craigslist$category[-s]) valid$document <- NULL # suppression du no du document (inutile) library(ranger) rg <- ranger(category ~ ., data=train, importance = "none", num.trees=300, mtry=20, replace=T, write.forest=T, seed=235) pred.rg <- predict(rg, dat = valid[,-ncol(valid)]) table(valid[,ncol(valid)], pred.rg$predictions) sum(pred.rg$predictions == valid[,ncol(valid)]) / nrow(valid) # autre forêt aléatoire (code absent du livre) # NB : temps de calcul très important library(randomForest) set.seed(235) rf <- randomForest(category ~ ., data=train, importance=FALSE, proximity=FALSE, ntree=300, mtry=20, replace=T, keep.forest=TRUE) pred.rf <- predict(rf, valid[,-ncol(valid)], type='prob') pred.rf[1:10,1:6] f <- function(i) { names(which.max(pred.rf[i,])) } pred.rdf <- Vectorize(f)(1:nrow(valid)) sum(pred.rdf == valid[,ncol(valid)]) / nrow(valid) table(valid[,ncol(valid)], pred.rdf) # boosting library(gbm) set.seed(235) boost <- gbm(category ~ ., data=train, n.trees = 500, distribution='multinomial', shrinkage = 0.2, n.minobsinnode = 5, verbose=T) pred.boost <- predict(boost, valid[,-ncol(valid)], n.trees=100, type="response") dim(pred.boost) f <- function(i) { names(which.max(pred.boost[i,,1])) } pred.gbm <- Vectorize(f)(1:nrow(valid)) sum(pred.gbm == valid[,ncol(valid)]) / nrow(valid) table(valid[,ncol(valid)], pred.gbm) # gradient boosting (code absent du livre) library(xgboost) library(Matrix) # transformation data frames en matrices sparses train.mx <- sparse.model.matrix(category ~ ., train) valid.mx <- sparse.model.matrix(category ~ ., valid) # transformation matrices sparses en matrices Xgb # à noter que la variable à expliquer doit être dans [0,1] (binary:logistic) ou dans [0,num_class] dtrain <- xgb.DMatrix(train.mx, label=as.numeric(train$category)-1) dvalid <- xgb.DMatrix(valid.mx, label=as.numeric(valid$category)-1) set.seed(235) system.time(train.gdbt <- xgb.train(params=list(objective="multi:softmax", num_class=6, eval_metric="mlogloss", eta=0.2, max_depth=5, colsample_bytree=0.5), data=dtrain, nrounds=1000, early_stopping_rounds = 10, watchlist=list(eval=dvalid))) pred.gbm <- predict(train.gdbt, newdata=dvalid) sum(pred.gbm == as.numeric(valid$category)-1) / nrow(valid) table(pred.gbm, as.numeric(valid$category)-1) # Section 10.4.3 : Pondération TF-IDF # --------------------------------------------------------------------------------------------------------- txt <- c("Au commencement était le Verbe, et le Verbe était en Dieu, et le Verbe était Dieu.", "Il était au commencement en Dieu.") library(quanteda) dtm <- dfm(txt) as.matrix(dtm) dtm <- dfm_tfidf(dtm) dtm <- dfm_tfidf(dtm, scheme_tf = "prop", scheme_df = "inverse") as.matrix(dtm) # Section 10.4.4 : Analyse sémantique latente (code absent du livre) # --------------------------------------------------------------------------------------------------------- (m <- matrix(c(1,2,0,2,1,3),2,3, byrow = TRUE)) (svd <- svd(m)) # matrice complète svd$u %*% diag(svd$d) %*% t(svd$v) # approximation de rang 1 diag(c(svd$d[1],0)) (a <- svd$u %*% diag(c(svd$d[1],0)) %*% t(svd$v)) a[2,]/a[1,] # LSA sur "La Recherche" de Proust library(proustr) proust <- as.data.frame(proust_books()) # analyses avec quanteda library(quanteda) corpus <- corpus(proust$text) # lemmatisation library(data.table) lemme <- fread("D:/Data Mining/Logiciel R/Programmes R/Text Mining/Lemmatisation/lemmes-fr.txt", sep="\t", header=F) setkey(lemme,V1) lemmatisation <- function(y) { unlist(lapply(lapply(strsplit(y,' '), function(x) as.character(ifelse(!is.na(lemme[x]$V2)&(lemme[x]$V3 %in% c("NOM", "ADJ")), lemme[x]$V2, ""))),function(x) paste(x,collapse =" "))) } texts(corpus) <- lemmatisation(texts(corpus)) # création matrice « documents x termes » dtm <- dfm(corpus, remove = unlist(proust_stopwords()), remove_punct = TRUE, stem = FALSE, tolower = TRUE) dtm[1:5,1:6] dtm # LSA avec package quanteda (plus rapide) lsa <- textmodel_lsa(dtm, nd=4) length(lsa$sk) # valeurs singulières (obtenues par SVD) dim(lsa$docs) # matrice des documents - distribution des thématiques par documents dim(lsa$features) # matrice des termes - distribution des mots par thématiques head(lsa$sk) # valeurs singulières (obtenues par SVD) head(lsa$docs) # matrice des documents - distribution des thématiques par documents head(lsa$features) # matrice des termes - distribution des mots par thématiques lsa$matrix_low_rank[1:5,1:6] # matrice de faible rang # calcul "à la main" de la matrice de faible rang lr <- lsa$docs %*% diag(lsa$sk) %*% t(lsa$features) lr[1:5,1:6] # 20 premiers termes de chaque thématique toptopics <- function(x){ top <- order(lsa$features[,x],decreasing=T) row.names(lsa$features)[head(top,20)] } toptopics(4) sapply(1:4,toptopics) # 1ère thématique de chaque document table(max.col(lsa$docs)) # matrice réduite dtm <- dfm_trim(dtm, min_termfreq = 300, min_docfreq = 30) dtm lsa <- textmodel_lsa(dtm, nd=4) library(ggplot2) library(ggrepel) # ACP sur matrice réduite pca <- prcomp(dtm) ggplot(as.data.frame(pca$rotation), aes(x=pca$rotation[,1], y=pca$rotation[,2])) + geom_point() + geom_text_repel(aes(label=rownames(pca$rotation)), hjust=1, vjust=1, size=3, segment.size = 0) # ACP sur matrice de faible rang pca <- prcomp(lsa$matrix_low_rank) ggplot(as.data.frame(pca$rotation), aes(x=pca$rotation[,1], y=pca$rotation[,2])) + geom_point() + geom_text_repel(aes(label=rownames(pca$rotation)), hjust=1, vjust=1, size=3, segment.size = 0) # AFC sur matrice réduite AFC <- textmodel_ca(dtm, nd = 2, sparse = TRUE) ggplot(as.data.frame(AFC$colcoord), aes(x=AFC$colcoord[,1], y=AFC$colcoord[,2])) + geom_point() + geom_text_repel(aes(label=rownames(AFC$colcoord)), hjust=1, vjust=1, size=3, segment.size = 0) # AFC sur matrice de faible rang (le graphique semble mieux séparer les termes que celui de l'ACP) AFC <- textmodel_ca(as.dfm(lsa$matrix_low_rank), nd = 2, sparse = TRUE) ggplot(as.data.frame(AFC$colcoord), aes(x=AFC$colcoord[,1], y=AFC$colcoord[,2])) + geom_point() + geom_text_repel(aes(label=rownames(AFC$colcoord)), hjust=1, vjust=1, size=3, segment.size = 0) # graphique des documents dk2 <- t(lsa$sk * t(lsa$docs)) head(dk2) plot(dk2[,1], y= dk2[,2], col="blue", pch="+", main="Documents Plot") text(dk2[,1], y= dk2[,2], labels=rownames(dk2), cex=.70) # graphique des termes (même graphique qu'une ACP) tk2 <- t(lsa$sk * t(lsa$features)) head(tk2) plot(tk2[,1], y= tk2[,2], col="red", cex=1, main="Terms Plot") text(tk2[,1], y= tk2[,2], labels=rownames(tk2) , cex=.80) # Section 10.4.5 : Fonction de densité de la loi de Dirichlet # --------------------------------------------------------------------------------------------------------- library(gtools) x1 <- seq(0.01,1, by=0.01) x2 <- seq(0.01,1, by=0.01) # Paramètre de concentration alpha <- c(0.1,0.1,0.1) alpha <- c(10,10,10) z <- matrix(0, nrow=length(x1), ncol=length(x2)) f <- function(x,y) { if(x+y < 1) {a <- c(x, y, 1-x-y) ; diric <- ddirichlet(a,alpha)} else {diric <- NA} return(diric) } z <- outer(x1,x2,Vectorize(f)) persp(x1,x2,z, xlab="alpha = (0.1,0.1,0.1)") persp(x1,x2,z, xlab="alpha = (10,10,10)") # Section 10.4.5 : LDA sur "La Recherche" de Proust # --------------------------------------------------------------------------------------------------------- library(proustr) proust <- as.data.frame(proust_books()) # analyses avec quanteda library(quanteda) length(unlist(tokens(proust[[1]]))) length(unlist(tokens(proust[[1]], remove_punct = TRUE))) table(proust$book, proust$year) corpus <- corpus(proust$text) # lemmatisation library(data.table) lemme <- fread("D:/Data Mining/Logiciel R/Programmes R/Text Mining/Lemmatisation/lemmes-fr.txt", sep="\t", header=F) setkey(lemme,V1) head(lemme);tail(lemme) lemmatisation <- function(y) { unlist(lapply(lapply(strsplit(y,' '), function(x) as.character(ifelse(!is.na(lemme[x]$V2)&(lemme[x]$V3 %in% c("NOM", "ADJ")), lemme[x]$V2, ""))),function(x) paste(x,collapse =" "))) } texts(corpus) <- lemmatisation(texts(corpus)) # création matrice « documents x termes » dtm <- dfm(corpus, remove = unlist(proust_stopwords()), remove_punct = TRUE, stem = FALSE, tolower = TRUE) dtm dtm[1:20,1:20] # LDA # ne peut être effectuée sur matrice DTM tm ou quanteda avec des lignes nulles (documents sans termes) # => nécessite conversion en "topicmodels" qui ne peut être effectuée que par quanteda ou suppression de ces lignes library(topicmodels) #lda <- LDA(dtm, k = 4, control = list(seed = 235)) # KO ! # suppression des lignes correspondant à des documents sans terme dtm <- dtm[rowSums(as.matrix(dtm)) > 0,] dtm <- convert(dtm, to = "topicmodels") # solution nécessaire pour calcul de perplexité # recherche de thématiques system.time(lda <- LDA(dtm, k = 4, control = list(seed = 235))) system.time(lda <- LDA(dtm, k = 4, control = list(seed = 235, alpha = 10))) system.time(lda <- LDA(dtm, k = 16, control = list(seed = 235), method = "Gibbs")) system.time(lda <- LDA(dtm, k = 4, control = list(seed = 235, burnin = 1000, iter = 1000, keep = 50), method = "Gibbs")) system.time(lda <- CTM(dtm, k = 4, control = list(seed = 235))) lda@control@alpha # paramètre de concentration logLik(lda) perplexity(lda, newdata=dtm) terms(lda,20) posterior(lda)$terms[1:4,1:10] t(posterior(lda)$terms[,1:10]) colnames(posterior(lda)$terms)[apply(posterior(lda)$terms,1,which.max)] # mot le plus probable de chaque topic table(topics(lda,1)) distHellinger(posterior(lda)$terms) # distance entre topics # recherche du nb optimal de thématiques ou de la meilleure méthode models <- list(VEM = LDA(dtm, k = 4, control = list(seed = 235), method = "VEM"), Gibbs = LDA(dtm, k = 4, control = list(seed = 235), method = "Gibbs"), CTM = CTM(dtm, k = 4, control = list(seed = 235))) head(posterior(models[[2]])$topics) # probas des thématiques dans chaque document sapply(models, function(x) mean(apply(posterior(x)$topics, 1, function(z) - sum(z * log(z))))) sapply(models[1:2], slot, "alpha") # paramètre de concentration (models.logLik <- as.matrix(lapply(models, logLik))) (models.perplexity <- as.matrix(lapply(models, function(x) perplexity(x, newdata=dtm)))) # on peut fixer la valeur du paramètre de concentration : lda <- LDA(dtm, k = 4, control = list(alpha = 1, seed = 235), method = "Gibbs") # recherche du nombre de thématiques optimal (code absent du livre) #install.packages("ldatuning") library(ldatuning) tunes <- FindTopicsNumber(dtm, topics = c(2:40), metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010"), method = "Gibbs", control = list(seed = 235), mc.cores = 4L, verbose = TRUE) tunes FindTopicsNumber_plot(tunes) # nombres de thématiques différents models <- lapply(seq(2, 6, by = 1), function(d){LDA(dtm, k=d, control = list(seed = 235), method = "Gibbs")}) sapply(models, function(x) mean(apply(posterior(x)$topics, 1, function(z) - sum(z * log(z))))) sapply(models, function(x) {mean(apply(posterior(x)$topics, 1, function(z) - sum(z * log(z))))/-log(1/x@k)}) # calculs d'entropie x <- c(0.5,0.5) x <- c(1/3,1/3,1/3) x <- c(0.1,0.2,0.7) x <- c(0.25,0.25,0.25,0.25) x <- c(0.2,0.2,0.3,0.3) - sum(x * log(x)) # Section 10.4.6 : Analyses de fréquences d’apparitions de termes # --------------------------------------------------------------------------------------------------------- # fréquence des termes (figure 10.5) library(proustr) proust <- as.data.frame(proust_books()) library(quanteda) corpus <- corpus(proust$text) corpus <- corpus(proust) corpus <- corpus_subset(corpus, book %in% c("Du côté de chez Swann")) dtm <- dfm(corpus, remove = unlist(proust_stopwords()), remove_punct = TRUE, stem = FALSE, tolower = TRUE, groups = "volume") keyness <- textstat_keyness(dtm, target = "Deuxième partie — Un amour de Swan") textplot_keyness(keyness) textplot_keyness(keyness, labelsize = 5, show_legend = F) kwic(corpus, "grand'mère") kwic(corpus, "prestige") # Section 10.4.7 : Word2Vec avec le package rword2vec # --------------------------------------------------------------------------------------------------------- library(devtools) install_github("mukul13/rword2vec") library(rword2vec) # apprentissage du modèle word2vec pour obtenir un vecteur de mots # remarque : text8 est le fichier d'apprentissage par défaut, constitué à partir de Wikipedia (anglais), # dont tous les mots sont séparés par un espace sans saut de ligne entre les articles # text8 est issu du corpus enwik9 (le premier milliard de caractères de Wikipédia) # http://mattmahoney.net/dc/text8.zip # http://mattmahoney.net/dc/enwik9.zip # http://mattmahoney.net/dc/textdata.html model <- word2vec(train_file = "D:/Data/NLP/text8", output_file = "text8.bin", binary=1) model <- word2vec(train_file = "D:/Data/NLP/enwik9.txt", output_file = "enwik9.bin", binary=1) # corpus français (http://fauconnier.github.io/) vecteur <- "D:/Data/NLP/frWac_non_lem_no_postag_no_phrase_500_skip_cut200.bin" #vecteur <- "D:/Data/NLP/frWiki_no_lem_no_postag_no_phrase_1000_skip_cut100.bin" #vecteur <- "D:/Data/NLP/frWiki_no_lem_no_postag_no_phrase_1000_cbow_cut100.bin" # corpus anglais (http://mattmahoney.net/dc/textdata.html) vecteur <- "D:/Data/NLP/enwik9.bin" # conversion du fichier binaire en fichier texte et affichage du vecteur de mots bin_to_txt(vecteur, "vector.txt") bin_to_txt("text8.bin", "vector.txt") w2v <- as.data.frame(read.table("vector.txt", skip=1)) w2v[1:10,] # distance (dist <- distance(file_name = "text8.bin", search_word = "king", num = 6)) (dist <- distance(file_name = vecteur, search_word = "roi", num = 20)) dist$word <- iconv(dist$word, from="UTF-8", to="latin1") dist # analogie (ana <- word_analogy(file_name = vecteur, search_words = "paris france berlin", num = 10)) (ana <- word_analogy(file_name = vecteur, search_words = "paris france rome", num = 10)) (ana <- word_analogy(file_name = vecteur, search_words = "homme femme roi", num = 10)) (ana <- word_analogy(file_name = vecteur, search_words = "man woman king", num = 6)) (ana <- word_analogy(file_name = "text8.bin", search_words = "man woman uncle", num = 6)) (ana <- word_analogy(file_name = vecteur, search_words = "beethoven mozart haydn", num = 10)) (ana <- word_analogy(file_name = vecteur, search_words = "san francisco los angeles", num = 10)) # Section 10.4.7 : Word2Vec avec le package wordVectors # --------------------------------------------------------------------------------------------------------- library(devtools) install_github("bmschmidt/wordVectors") library(wordVectors) library(magrittr) vecteur <- "D:/Data/NLP/text8.bin" vecteur <- "D:/Data/NLP/enwik9.bin" vecteur <- "D:/Data/NLP/frWac_non_lem_no_postag_no_phrase_500_skip_cut200.bin" # lecture d'un modèle existant model <- read.vectors(vecteur) # création d'un modèle à partir d'un corpus model <- train_word2vec("D:/Data/NLP/craigswords.csv", "craigs2.bin", vectors=100, threads=4, window=5, iter=10, cbow=0, min_count=5, negative_samples=0) #if (!file.exists("corpus_vectors.bin")) {model = train_word2vec("corpus.txt","corpus_vectors.bin",vectors=200,threads=4,window=12,iter=10,negative_samples=0)} else model = read.vectors("corpus_vectors.bin") # positions de mots model[["roi"]]@.Data # coordonnées d'un mot (model[["king"]]+model[["queen"]])/2 # moyenne des coordonnées de deux mots model %>% closest_to("roi") (ana <- closest_to(model, model[["king"]])) (ana <- closest_to(model, model[[c("king","queen")]])) model %>% closest_to(model[[c("king","man","woman")]], 6) model %>% closest_to(model[["king"]]-model[["man"]]+model[["woman"]], 6) model %>% closest_to(model[["roi"]]-model[["homme"]]+model[["femme"]], 6) model %>% closest_to(~"king" - "man" + "woman") model %>% closest_to(~"roi" - "homme" + "femme", 5) (ana <- closest_to(model, model[[c("paris","londres","berlin")]], 50)) (ana <- closest_to(model, model[[c("roi","reine")]], 50)) groupe <- model[[ana$word, average=F]] plot(groupe, method="pca") # Section 10.4.7 : Word2Vec : Exemple de détection des spams # --------------------------------------------------------------------------------------------------------- library(h2o) h2o.init() # lecture d'un fichier de spam récupéré sur Kaggle https://www.kaggle.com/uciml/sms-spam-collection-dataset classes <- c("Enum", "String") textes <- h2o.importFile(path="D:/Data/NLP/data_spam.csv", header = TRUE, col.types = classes) head(textes) nrow(textes) # Tokenisation du texte en mots words <- h2o.tokenize(textes$v2, "\\\\W+") dim(words) head(as.data.frame(words), 20) # Construction du modèle Word2Vec system.time(w2v.model <- h2o.word2vec(words, min_word_freq = 5, vec_size = 100, window_size = 5, epochs = 10)) # Synonymes d'un mot h2o.findSynonyms(w2v.model, "free", count = 10) # Calcul des vecteurs Word2Vec plot.vecs <- h2o.transform(w2v.model, words, aggregate = "average") # Création des échantillons d'apprentissage et test data <- h2o.cbind(textes["v1"], plot.vecs) data.split <- h2o.splitFrame(data, ratios = 0.8, seed = 235) head(data.split) h2o.table(data.split[[2]]$v1) # Modèle de gradient boosting gbm.model <- h2o.gbm(x = names(plot.vecs), y = "v1", training_frame = data.split[[1]], validation_frame = data.split[[2]], seed = 235) gbm.model <- h2o.gbm(x = names(plot.vecs), y = "v1", training_frame = data.split[[1]], validation_frame = data.split[[2]], ntrees = 100, learn_rate = 0.2, max_depth = 5, seed = 235) summary(gbm.model) h2o.varimp(gbm.model) # matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(gbm.model, data.split[[2]]) # taux d'erreur issu de la matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(gbm.model, data.split[[2]])$Error[3] # Modèle forêts aléatoires rf.model <- h2o.randomForest(x = names(plot.vecs), y = "v1", training_frame = data.split[[1]], validation_frame = data.split[[2]], ntrees = 100, max_depth = 20, mtries = 10, seed=235) h2o.confusionMatrix(rf.model, data.split[[2]]) # Modèle réseau de neurones dl.model <- h2o.deeplearning(x = names(plot.vecs), y = "v1", training_frame = data.split[[1]], validation_frame = data.split[[2]], hidden = c(10,10), epochs = 20, activation = "Tanh", l1 = 1e-6, seed=235) h2o.confusionMatrix(dl.model, data.split[[2]]) # Section 10.4.7 : Word2Vec : Exemple de détection du type d'offre d'emploi # --------------------------------------------------------------------------------------------------------- library(h2o) h2o.init() jobs.path = "https://raw.githubusercontent.com/h2oai/sparkling-water/rel-1.6/examples/smalldata/craigslistJobTitles.csv" jobs <- h2o.importFile(jobs.path, destination_frame = "jobs", col.names = c("category", "jobtitle"), col.types = c("Enum", "String"), header = TRUE) head(jobs,10) nrow(jobs) h2o.table(jobs$category) # découpage des titres d'annonces en mots - version courte (sans preprocessing) words <- h2o.tokenize(jobs$jobtitle, "\\\\W+") # découpage des titres d'annonces en mots - version longue library(tm) # pour liste des stop-words tokenize <- function(sentences) { tokenized <- h2o.tokenize(sentences, "\\\\W+") # convert to lower case tokenized.lower <- h2o.tolower(tokenized) # remove short words (less than 2 characters) tokenized.lengths <- h2o.nchar(tokenized.lower) tokenized.filtered <- tokenized.lower[is.na(tokenized.lengths) || tokenized.lengths >= 2,] # remove words that contain numbers tokenized.words <- tokenized.filtered[h2o.grep("[0-9]", tokenized.filtered, invert = TRUE, output.logical = TRUE),] # remove stop words #tokenized.words[is.na(tokenized.words) || (! tokenized.words %in% stop_words),] tokenized.words[(! tokenized.words %in% stopwords("en")),] } # découpage des titres d'annonces en mots words <- tokenize(jobs$jobtitle) nrow(words) head(words,10) # modèle word2vec w2v.model <- h2o.word2vec(words, min_word_freq = 5, vec_size = 100, window_size = 5, epochs = 10) # vérification (synonymes de "teacher") print(h2o.findSynonyms(w2v.model, "teacher", count = 5)) # calcul des vecteurs Word2Vec jobs.vecs <- h2o.transform(w2v.model, words, aggregate = "average") # création des échantillons d'apprentissage et test data <- h2o.cbind(jobs$category, jobs.vecs) data <- data[!is.na(data$C1),] data.split <- h2o.splitFrame(data, ratios = 0.8, seed = 235) head(data.split) h2o.table(data.split[[2]]$category) # modèle de gradient boosting gbm.model <- h2o.gbm(x = names(jobs.vecs), y = "category", training_frame = data.split[[1]], validation_frame = data.split[[2]], distribution="multinomial", ntrees = 300, learn_rate = 0.1, max_depth = 6, seed = 235) summary(gbm.model) # taux d'erreur issu de la matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(gbm.model, data.split[[2]])$Error[7] # matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(gbm.model, data.split[[2]]) # Modèle forêts aléatoires rf.model <- h2o.randomForest(x = names(jobs.vecs), y = "category", training_frame = data.split[[1]], validation_frame = data.split[[2]], ntrees = 200, max_depth = 30, mtries = 10, seed=235) # taux d'erreur issu de la matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(rf.model, data.split[[2]])$Error[7] # matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(rf.model, data.split[[2]]) # Modèle réseau de neurones (code absent du livre) dl.model <- h2o.deeplearning(x = names(jobs.vecs), y = "category", training_frame = data.split[[1]], validation_frame = data.split[[2]], hidden = c(10,10), epochs = 10, activation = "Rectifier", max_w2 = 10, seed=235) # taux d'erreur issu de la matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(dl.model, data.split[[2]])$Error[7] # matrice de confusion sur l'échantillon de validation h2o.confusionMatrix(dl.model, data.split[[2]]) # représentation des mots # création d'un data frame contenant chaque mot (répété autant de fois qu'il l'est dans le corpus) et ses vecteurs mots.vecs <- h2o.transform(w2v.model, words, aggregate = "none") words$mot <- words$C1 # renommage variable car C1 est aussi le nom du 1er vecteur Word2Vec words$C1 <- NULL mots.vecs <- h2o.cbind(words, mots.vecs) mots.vecs <- mots.vecs[!is.na(mots.vecs$C1),] nrow(mots.vecs) # sélection des mots-vecteurs les plus fréquents motsvectsdf <- as.data.frame(mots.vecs) motsvectsdf <- unique(motsvectsdf) # suppression des mots en double nrow(motsvectsdf) freq <- as.data.frame(words$mot) freq <- table(freq$mot) freq <- sort(freq, decreasing=TRUE) select <- names(head(freq, 50)) select # sélection des mots-vecteurs les plus fréquents mots.vect <- motsvectsdf[motsvectsdf$mot %in% select, c("mot","C1","C2")] library(ggplot2) ggplot(mots.vect, aes(x=C1, y=C2)) + geom_point() + geom_text(aes(label=mot), hjust=1, vjust=1, size=5) library(ggrepel) ggplot(mots.vect, aes(x=C1, y=C2)) + geom_point() + geom_text_repel(aes(label=mot), hjust=1, vjust=1, size=8, segment.size = 0) # plan factoriel pca <- prcomp(motsvectsdf[,-1]) pcamots <- cbind(motsvectsdf$mot, pca$x[,1:2]) colnames(pcamots)[1] <- "mot" motsvectspca <- as.data.frame(pcamots) motsvectspca$PC1 <- as.numeric(as.character(motsvectspca$PC1)) motsvectspca$PC2 <- as.numeric(as.character(motsvectspca$PC2)) mots.vect <- motsvectspca[motsvectspca$mot %in% select, ] ggplot(as.data.frame(mots.vect), aes(x=PC1, y=PC2)) + geom_point() + geom_text(aes(label=mot), hjust=1, vjust=1, size=5) ggplot(as.data.frame(mots.vect), aes(x=PC1, y=PC2)) + geom_point() + geom_text_repel(aes(label=mot), hjust=1, vjust=1, size=8, segment.size = 0) # fermeture du cluster H2O h2o.shutdown(prompt=FALSE) # Section 10.4.8 : Représentation GloVe (code absent du livre) # --------------------------------------------------------------------------------------------------------- install.packages('text2vec') library(text2vec) detach("package:h2o") browseVignettes(package = "text2vec") # récupération d'un corpus en local corpus <- "D:/Data/NLP/text8" wiki <- readLines(corpus, n = 1, warn = FALSE) # création du vocabulaire (de mots = unigrammes) tokens <- space_tokenizer(wiki) # + rapide tokens <- word_tokenizer(wiki) # + général mais plus lent vocab <- create_vocabulary(itoken(tokens)) # suppression des mots rares (apparaissant moins de 5 fois) vocab <- prune_vocabulary(vocab, term_count_min = 5, doc_proportion_max = 0.5) # construction de la matrice de co-occurrence # à partir du vocabulaire filtré vectorizer <- vocab_vectorizer(vocab) # fenêtre de contexte de largeur 5 tcm <- create_tcm(it, vectorizer, skip_grams_window = 5) tcm[1:20,1:20] # apprentissage du modèle glove = GlobalVectors$new(rank = 50, x_max = 10) glove = GloVe$new(rank = 50, x_max = 10) wv_main = glove$fit_transform(tcm, n_iter = 10, convergence_tol = 0.01, n_threads = 8) dim(wv_main) library(Matrix) rankMatrix(wv_main) wv_main[1:30,1:15] wv_context = glove$components dim(wv_context) rankMatrix(wv_context) # get average of main and context vectors as proposed in GloVe paper word_vectors = wv_main + t(wv_context) dim(word_vectors) word_vectors[1:30,1:15] berlin <- word_vectors["paris", , drop = FALSE] - word_vectors["france", , drop = FALSE] + word_vectors["germany", , drop = FALSE] cos_sim = sim2(x = word_vectors, y = berlin, method = "cosine", norm = "l2") head(sort(cos_sim[,1], decreasing = TRUE), 5) # Section 10.4.9 : fastText : Exemple de détection du type d'offre d'emploi # --------------------------------------------------------------------------------------------------------- #install.packages('fastTextR') library(fastTextR) # exemple Kaggle de spams spams <- readLines("D:/Data/NLP/base_spam.csv") head(spams,10) spams <- normalize(spams) set.seed(123) s <- sort(sample(length(spams), length(spams)*0.8, replace=F)) table(gsub("\\D", "", substr(spams[s],1,10))) writeLines(spams[s], con = "D:/Data/NLP/spam.train") writeLines(spams[-s], con = "D:/Data/NLP/spam.test") cntrl <- ft.control(word_vec_size = 100, window_size = 5, learning_rate = 0.05, min_count = 5, epoch = 10, nthreads = 4) cntrl <- ft.control(word_vec_size = 100, window_size = 5, learn_update = 100, learning_rate = 0.05, min_count = 5, min_ngram = 2, max_ngram = 6, epoch = 10, nthreads = 2) system.time(model <- fasttext(input="D:/Data/NLP/spam.train", method="supervised", control = cntrl)) model # application du modèle test.pred <- predict(model, newdata_file="D:/Data/NLP/spam.test", k = 1, prob = TRUE) (confusion_matrix <- table(gsub("\\D", "", substr(spams[-s],1,10)), gsub("\\D", "", test.pred$label))) 1-sum(diag(confusion_matrix)) / sum(confusion_matrix) # exemple de la CraigsList library(h2o) localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, max_mem_size = '4g', ice_root="D:/Data/h2o", nthreads=2) jobs.path = "https://raw.githubusercontent.com/h2oai/sparkling-water/rel-1.6/examples/smalldata/craigslistJobTitles.csv" jobs <- h2o.importFile(jobs.path, destination_frame = "jobs", col.names = c("category", "jobtitle"), col.types = c("Enum", "String"), header = TRUE) jobs <- as.data.frame(jobs) jobs$category <- paste0("__label__", jobs$category) head(jobs,20) # prétraitement jobs$jobtitle <- gsub("[^[:print:]]", "", jobs$jobtitle) # suppression des caractères non imprimables #jobs$jobtitle <- gsub("\\<\\w{1,2}\\>", "", jobs$jobtitle) # suppression des mots de 1 ou 2 caractères jobs$jobtitle <- gsub("\\<\\w{1}\\>", "", jobs$jobtitle) # suppression des mots de 1 caractère jobs$jobtitle <- tolower(jobs$jobtitle) # passage en bas de casse jobs$jobtitle <- gsub("[[:punct:]]", " ", jobs$jobtitle) # remplacement de la ponctuation par un espace jobs$jobtitle <- gsub("\\<[a-z]*[0-9]+[a-z]*\\>", "", jobs$jobtitle) # suppression des mots contenant un chiffre #jobs$jobtitle <- gsub("[[:digit:]]", "", jobs$jobtitle) # suppression des chiffres # suppression des mots-outils library(tm) stopw <- paste0("\\b(", paste0(stopwords("en"), collapse="|"), ")\\b") #stopw <- paste0("\\b(", paste0(list.words, collapse="|"), ")\\b") jobs$jobtitle <- gsub(stopw, "", jobs$jobtitle) # suppression des mots-outils jobs$jobtitle <- gsub("[[:space:]]{2,}", " ", jobs$jobtitle) # suppression des espaces multiples # sélection des mots les plus fréquents (à faire avant de modifier "jobs" 10 lignes plus bas par readLines) head(jobs, 20) mon_tokenizer <- function(x) unlist(strsplit(as.character(x), "[[:space:]]+")) freq <- mon_tokenizer(jobs$jobtitle) freq <- table(freq) freq <- sort(freq, decreasing=TRUE) select <- names(head(freq, 50)) select # enregistrement du fichier pour utilisation par fastText write.table(jobs[-1,], file="D:/Data/NLP/craigsjobs.csv", row.names = F, col.names = F) jobs <- readLines("D:/Data/NLP/craigsjobs.csv") jobs <- normalize(jobs) # modélisation set.seed(123) s <- sample(length(jobs), length(jobs)*0.8, replace=F) writeLines(jobs[s], con = "D:/Data/NLP/craigslist.train") writeLines(jobs[-s], con = "D:/Data/NLP/craigslist.test") cntrl <- ft.control(word_vec_size = 20, window_size = 5, learn_update = 100, learning_rate = 0.1, min_count = 5, min_ngram = 2, max_ngram = 5, epoch = 10, nthreads = 4) system.time(model <- fasttext(input="D:/Data/NLP/craigslist.train", method="supervised", control = cntrl)) model # application du modèle test.pred <- predict.supervised_model(model, newdata_file="D:/Data/NLP/craigslist.test", k = 1, prob = TRUE) table(test.pred$label) (confusion_matrix <- table(gsub("__label__", "", substr(jobs[-s],1,21)), gsub("__label__", "", test.pred$label))) 1-sum(diag(confusion_matrix)) / sum(confusion_matrix) # représentation des mots de la CraigsList cntrl <- ft.control(learning_rate = 0.1, word_vec_size = 20, epoch = 10, nthreads = 4) model <- fasttext("D:/Data/NLP/craigslist.train", "skipgram", cntrl) # vecteurs des mots de la Craigslist vecteurs <- get_word_vectors(model, get_words(model)) dim(vecteurs) # mise en forme des mots-vecteurs motsvectsdf <- as.data.frame(vecteurs) motsvectsdf$mot <- row.names(motsvectsdf) # sélection des mots-vecteurs les plus fréquents mots.vect <- motsvectsdf[motsvectsdf$mot %in% select, ] library(ggplot2) ggplot(mots.vect, aes(x=V1, y=V2)) + geom_point() + geom_text(aes(label=mot), hjust=1, vjust=1, size=5) library(ggrepel) ggplot(mots.vect, aes(x=V1, y=V2)) + geom_point() + geom_text_repel(aes(label=mot), hjust=1, vjust=1, size=8, segment.size = 0) # représentation dans le plan factoriel pca <- prcomp(motsvectsdf[,-ncol(motsvectsdf)]) pcamots <- cbind(motsvectsdf$mot,pca$x[,1:2]) colnames(pcamots)[1] <- "mot" motsvectspca <- as.data.frame(pcamots) motsvectspca$PC1 <- as.numeric(as.character(motsvectspca$PC1)) motsvectspca$PC2 <- as.numeric(as.character(motsvectspca$PC2)) mots.vect <- motsvectspca[motsvectspca$mot %in% select, ] ggplot(as.data.frame(mots.vect), aes(x=PC1, y=PC2)) + geom_point() + geom_text(aes(label=mot), hjust=1, vjust=1, size=5) ggplot(as.data.frame(mots.vect), aes(x=PC1, y=PC2)) + geom_point() + geom_text_repel(aes(label=mot), hjust=1, vjust=1, size=8, segment.size = 0) # Section 10.5.1 : génération de texte # --------------------------------------------------------------------------------------------------------- library(keras) library(readr) library(stringr) library(tokenizers) # texte source # http://www.gutenberg.org/ebooks/29843 path <- "D:/Data/NLP/contemplations_autrefois.txt" # nombre de mots dans le texte source text <- read_lines(path) %>% str_to_lower() %>% str_c(collapse = "\n") library(quanteda) corpus <- corpus(text) corpus <- corpus(texts(corpus, groups = rep(1, ndoc(corpus)))) ntoken(corpus) ntype(corpus) ntoken(corpus) / ntype(corpus) # paramètres maxlen <- 64 step <- 1 # Préparation des données -------------------------------------------------------- # On lit le texte, on passe toutes les capitales en bas de casse, on conserve la ponctuation et les espaces, et on découpe (tokenise) le texte en caractères text <- read_lines(path) %>% str_to_lower() %>% str_c(collapse = "\n") %>% tokenize_characters(strip_non_alphanum = FALSE, simplify = TRUE) print(sprintf("corpus length: %d", length(text))) chars <- text %>% unique() %>% sort() print(sprintf("total chars: %d", length(chars))) # On découpe le texte en séquences recouvrantes de maxlen caractères library(purrr) dataset <- map(seq(1, length(text) - maxlen - 1, by = step), ~list(sentence = text[.x:(.x + maxlen - 1)], next_char = text[.x + maxlen])) dataset <- transpose(dataset) # Codage des caractères en tenseurs binaires x <- array(0L, dim = c(length(dataset$sentence), maxlen, length(chars))) y <- array(0L, dim = c(length(dataset$sentence), length(chars))) for(i in 1:length(dataset$sentence)){ x[i,,] <- sapply(chars, function(x){as.integer(x==dataset$sentence[[i]])}) y[i,] <- as.integer(chars==dataset$next_char[[i]]) } # Définition du modèle LSTM prédisant le prochain caractère # (la sortie softmax a autant d'unités qu'il y a de caractères différents à prédire) model <- keras_model_sequential() model %>% layer_lstm(units = 512, return_sequences=TRUE, input_shape = c(maxlen, length(chars))) %>% layer_dropout(0.25, seed=235) %>% layer_lstm(units = 512, return_sequences=TRUE) %>% layer_dropout(0.25, seed=235) %>% layer_lstm(units = 512) %>% layer_dropout(0.25, seed=235) %>% layer_dense(length(chars)) %>% layer_activation("softmax") optimizer <- optimizer_rmsprop(lr = 0.001) model %>% compile(loss = "categorical_crossentropy", optimizer = optimizer) model %>% compile(loss = "categorical_crossentropy", optimizer = optimizer, metrics = c("accuracy")) # Tirage du prochain caractère selon les prédictions du modèle LSTM sample_next_char <- function(preds, temperature = 1){ preds <- log(preds)/temperature exp_preds <- exp(preds) preds <- exp_preds/sum(exp(preds)) which.max(t(rmultinom(1, 1, preds))) } # Génération d'un texte de 500 caractères on_epoch_end <- function(epoch, logs) { start_index <- sample(1:(length(text) - maxlen), size = 1) #cat(sprintf("epoch: %02d ---------------\n\n", epoch)) # test de plusieurs températures (entropies) for(temperature in c(0.1, 0.3, 0.5)){ cat(sprintf("température : %f ---------------\n", temperature)) # tirage au hasard d'un échantillon de texte de longueur maxlen sentence <- text[start_index:(start_index + maxlen - 1)] generated <- "" cat("graine de texte : ", sentence, "\n\n") # génération de 500 caractères : on code l'échantillon d'apprentissage, on applique le modèle LSTM, # on prédit le prochain caractère (en tenant compte de la température), on l'ajoute au texte généré # et on continue en décalant l'échantillon d'un caractère for(i in 1:500){ s <- sapply(chars, function(s){as.integer(s==sentence)}) s <- array_reshape(s, c(1, dim(s))) # prédiction du prochain caractère selon le modèle preds <- predict(model, s) next_index <- sample_next_char(preds, temperature) next_char <- chars[next_index] generated <- paste0(generated, next_char) # mise à jour de l'échantillon d'apprentissage sentence <- c(sentence[-1], next_char) } cat(generated) cat("\n\n") } } rm(dataset) gc() summary(model) generate_text <- callback_lambda(on_epoch_end = on_epoch_end) set.seed(123) system.time(model %>% fit(x, y, batch_size = 128, epochs = 20, callbacks = generate_text)) callbacks_list <- list(callback_early_stopping(monitor = "accuracy", patience = 1), callback_model_checkpoint(filepath = "hugo_lstm_model.h5", monitor = "val_loss", save_best_only = TRUE), callback_lambda(on_epoch_end = on_epoch_end)) set.seed(123) system.time(model %>% fit(x, y, batch_size = 128, epochs = 30, callbacks = callbacks_list)) # enregistrement du modèle save_model_hdf5(model, "hugo-512-512-512_lstm_model.h5") # modèle entier enregistré comme fichier binaire save_model_weights_hdf5(model, "weights_hugo-512-512-512_lstm_model.h5") # ensemble des poids du modèle enregistré comme fichier binaire # Sections 10.5.2 et 10.5.4 : LSTM et CNN pour classement de plaintes # --------------------------------------------------------------------------------------------------------- library(data.table) system.time(textes <- fread("D:/Data/NLP/Consumer_Complaints.csv", check.names=TRUE)) setnames(textes, old = "Consumer.complaint.narrative", new = "text") dim(textes) head(textes) textes <- textes[, .(Product,text)] textes <- textes[text!=""] table(textes$Product) # prétraitement textes$text <- gsub("[^[:print:]]", "", textes$text) # suppression des caractères non imprimables textes$text <- gsub("\\<\\w{1,2}\\>", "", textes$text) # suppression des mots de 1 ou 2 caractères textes$text <- gsub("\\<[a-z]*[x]{2,}[a-z]*\\>", "", textes$text) # suppression des mots contenant une chaîne de xx textes$text <- gsub("\\<[a-z]*[X]{2,}[a-z]*\\>", "", textes$text) # suppression des mots contenant une chaîne de XX textes$text <- gsub("[[:punct:]]", " ", textes$text) # suppression de la ponctuation (remplacement par un espace) textes$text <- gsub("[[:space:]]{2,}", " ", textes$text) # suppression des espaces multiples # exclusion des deux produits les plus rares textes <- textes[!(textes$Product%in%c("Other financial service","Virtual currency")),] dim(textes) head(textes) table(textes$Product) # préparation des données X <- textes$text Y <- textes$Product Y <- as.numeric(factor(Y)) set.seed(235) s <- sample(1:length(X),length(X)*0.7) x_train <- X[s] x_test <- X[-s] y_train <- Y[s] y_test <- Y[-s] maxwords <- 5000 # 10000 maxlength <- 200 # 500 (num_classes <- length(table(y_train))) # prédictions aléatoires pred <- sample(y_test) # précision mean(y_test == pred) rm(pred) # modèle LSTM library(keras) library(tokenizers) # échantillon d'apprentissage # découpage du texte en mots, dont le résultat est une liste dont chaque élément est un document du corpus #x_train <- tokenize_words(x_train, lowercase = TRUE, stopwords = stopwords::stopwords("en"), strip_punct = TRUE, strip_numeric = FALSE, simplify = FALSE) x_token <- tokenize_words(x_train, lowercase = TRUE, stopwords = NULL, strip_punct = TRUE, strip_numeric = FALSE, simplify = FALSE) head(x_token,1) mean(sapply(1:length(x_token),function(x) length(x_token[[x]]))) # nombre de mots moyen par message quantile((sapply(1:length(x_token),function(x) length(x_token[[x]]))),0.8) # nombre de mots moyen par message # transformation de chaque mot en un entier qui est le numéro d'un mot entre 1 et maxwords # deux mots différents peuvent avoir le même numéro si maxwords est inférieur au nombre total de mots différents x_encode <- lapply(x_token, function(x) text_one_hot(paste(x, collapse=" "), maxwords)) head(x_encode,1) #x_encode <- lapply(x_token, function(x) text_one_hot(x, maxwords, lower=FALSE, filters = "!\"#$%&()*+,-./:;<=>?@[\\]^_`{|}~\t\n")) # cette fonction transforme une liste de n suites d'entiers (ici les numéros de mots) en une matrice de dimension (n,maxlen) # les séquences plus courtes que maxlen sont complétées par des 0 au début de la séquence x_train <- pad_sequences(x_encode, maxlen = maxlength) dim(x_train) head(x_train,1) # échantillon de test x_token <- tokenize_words(x_test, lowercase = TRUE, stopwords = NULL, strip_punct = TRUE, strip_numeric = FALSE, simplify = FALSE) # transformation de chaque mot en un entier qui est le numéro d'un mot entre 1 et maxwords # deux mots différents peuvent avoir le même numéro si maxwords est inférieur au nombre total de mots différents x_encode <- lapply(x_token, function(x) text_one_hot(paste(x, collapse=" "), maxwords)) # cette fonction transforme une liste de n suites d'entiers (ici les numéros de mots) en une matrice de dimension (n,maxlen) # les séquences plus courtes que maxlen sont complétées par des 0 au début de la séquence x_test <- pad_sequences(x_encode, maxlen = maxlength) dim(x_test) y_train <- to_categorical(y_train-1, num_classes) y_test <- to_categorical(y_test-1, num_classes) # définition du modèle LSTM model <- keras_model_sequential() model %>% #layer_lstm(units = 128, input_shape = c(maxlen, 100)) %>% layer_embedding(input_dim = maxwords, output_dim = 100, input_length=maxlength) %>% #layer_lstm(units = 128, return_sequences=TRUE) %>% #layer_dropout(0.2, seed=235) %>% #layer_lstm(units = 128) %>% layer_lstm(units = 128, dropout=0.2, recurrent_dropout = 0.2) %>% #layer_dropout(0.2, seed=235) %>% layer_dense(units = 64, activation = "relu") %>% #layer_batch_normalization(axis=1) %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle GRU model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 100, input_length=maxlength) %>% layer_gru(units = 128, recurrent_dropout = 0.2) %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle GRU bidirectionnel model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 100, input_length=maxlength) %>% bidirectional(layer_gru(units = 128)) %>% layer_dropout(0.2, seed=235) %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle CNN-LSTM model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 100, input_length=maxlength) %>% layer_conv_1d(filters = 64, kernel_size = 7, activation = "relu") %>% layer_max_pooling_1d(pool_size = 5) %>% layer_conv_1d(filters = 64, kernel_size = 7, activation = "relu") %>% layer_lstm(units = 128) %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle CNN model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 100, input_length=maxlength) %>% layer_conv_1d(filters = 128, kernel_size = 7, activation = "relu") %>% layer_max_pooling_1d(pool_size = 5) %>% layer_conv_1d(filters = 128, kernel_size = 7, activation = "relu") %>% layer_global_max_pooling_1d() %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle CNN (code absent du livre) model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 100, input_length=maxlength) %>% layer_conv_1d(filters = 64, kernel_size = 7) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_max_pooling_1d(pool_size = 5) %>% layer_conv_1d(filters = 64, kernel_size = 7) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_global_max_pooling_1d() %>% # max pooling global ou aplatissement layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # compilation model %>% compile(loss = 'categorical_crossentropy', optimizer = 'rmsprop', metrics = c('accuracy')) # apprentissage #model %>% fit(x_train, y_train, batch_size = 64, epochs = 10) system.time(history <- fit(model, x_train, y_train, batch_size = 128, epochs = 10, shuffle = TRUE, validation_data = list(x_test, y_test))) # prédictions pred <- predict_classes(model, x_test) test.label <- max.col(y_test) - 1 table(test.label, pred) mean(test.label != pred) # historique de l'apprentissage plot(history) # Section 10.5.3 : Classement de textes à l’aide d’un modèle H2O # --------------------------------------------------------------------------------------------------------- library(h2o) h2o.init(max_mem_size = '6g', ice_root="D:/Data/h2o", nthreads=2) system.time(textes <- h2o.importFile(path = "D:/Data/NLP/Complaints.csv", header=TRUE, sep=",")) # nb de modalités à prédire nrow(h2o.table(textes$Product)) textes <- textes[!(textes$Product%in%c("Other financial service","Virtual currency")),] # Tokenisation du texte en mots words <- h2o.tokenize(textes$text, "\\\\W+") dim(words) # Construction du modèle Word2Vec system.time(w2v.model <- h2o.word2vec(words, min_word_freq = 50, vec_size = 100, window_size = 5, epochs = 10)) # Synonymes d'un mot h2o.findSynonyms(w2v.model, "complaint", count = 10) # Calcul des vecteurs Word2Vec system.time(plot.vecs <- h2o.transform(w2v.model, words, aggregate = "average")) dim(plot.vecs) # Création des échantillons d'apprentissage et test data <- h2o.cbind(textes["Product"], plot.vecs) data.split <- h2o.splitFrame(data, ratios = 0.7, seed = 235) # réseau de neurones # grid search dl_params <- list(activation = c("Tanh","Rectifier"), l1 = c(1e-4,1e-5,1e-6,0), l2 = c(1e-4,1e-5,1e-6,0), hidden = c(c(25,25),c(50,50),c(100,100),c(25,25,25),c(50,50,50))) system.time(grid_dl <- h2o.grid(algorithm="deeplearning", grid_id="dl_grid", epochs=20, training_frame=data.split[[1]], validation_frame=data.split[[2]], x=names(plot.vecs), y="Product", seed=123, hyper_params=dl_params)) grid_dl # modèles triés par erreur croissante en validation models <- h2o.getGrid("dl_grid", sort_by="err", decreasing=FALSE) # paramètres de la grille et leurs valeurs pour les meilleurs et pires modèles models models@summary_table # paramètres du meilleur modèle models@summary_table[1,] # paramètres avec en plus liste décroissante des modèles summary(models) # extraction du meilleur modèle (selon tri précédent) best_model <- h2o.getModel(models@model_ids[[1]]) best_model # gradient boosting # grid search gbm_params <- list(ntrees = c(100, 300), learn_rate = c(0.2,0.35,0.5), max_depth = c(5,10)) system.time(grid_gbm <- h2o.grid(algorithm="gbm", grid_id="dl_grid", training_frame=data.split[[1]], validation_frame=data.split[[2]], x=names(plot.vecs), y="Product", seed=123, hyper_params=gbm_params)) grid_gbm # modèles triés par erreur croissante en validation models <- h2o.getGrid("dl_grid", sort_by="err", decreasing=FALSE) # paramètres de la grille et leurs valeurs pour les meilleurs et pires modèles models models@summary_table # paramètres du meilleur modèle models@summary_table[1,] # paramètres avec en plus liste décroissante des modèles summary(models) # extraction du meilleur modèle (selon tri précédent) best_model <- h2o.getModel(models@model_ids[[1]]) best_model # Section 10.5.5 : LSTM pour détection de spams # --------------------------------------------------------------------------------------------------------- # lecture à partir d'un fichier lu par h2o et enregistré comme fichier à deux variables # en effet, la lecture directe du fichier Kaggle par read.table ou fread est impossible # car le séparateur est la "," mais qu'il peut y avoir plusieurs "," sur une ligne # et que seule la première virgule doit être considérée comme un séparateur library(data.table) spams <- fread("D:/Data/NLP/emails_spam.csv", head=F, sep=";", col.names = c("label","texte"), colClasses = c("factor","character")) dim(spams) str(spams) table(spams$label) # données de spams X <- spams$texte Y <- spams$label Y <- as.numeric(factor(Y)) # échantillons set.seed(235) s <- sample(1:length(X),length(X)*0.8) x_train <- X[s] x_test <- X[-s] y_train <- Y[s] y_test <- Y[-s] (num_classes <- length(table(y_train))) # nombre de mots distincts library(quanteda) length(unique(unlist(tokens(spams$texte, remove_punct = TRUE)))) maxwords <- 10000 # modèle LSTM library(keras) library(tokenizers) # échantillon d'apprentissage # découpage du texte en mots, dont le résultat est une liste dont chaque élément est un document du corpus, i.e. un mail x_token <- tokenize_words(x_train, lowercase = FALSE, stopwords = NULL, strip_punct = TRUE, strip_numeric = FALSE) # transformation de chaque mot en un entier qui est le numéro d'un mot entre 1 et maxwords # deux mots différents peuvent avoir le même numéro si maxwords est inférieur au nombre total de mots différents # mais ce codage s’appuie sur une fonction de hash optimisée x_encode <- lapply(x_token, function(x) text_one_hot(paste(x, collapse=" "), maxwords)) # cette fonction ne s’applique pas à chaque élément x mais à paste(x, collapse=" "), # car cette fonction s’applique à des chaînes de mots « aaa bbb ccc… » et non à des vecteurs (« aaa », « bbb », « ccc »…) # la fonction pad_sequences transforme une liste de n suites d'entiers (ici les numéros de mots) en une matrice de dimension (n,maxlen) # les séquences plus courtes que maxlen sont complétées par des 0 au début de la séquence # les séquences plus longues sont tronquées # on a fixé ce paramètre maxlength à 25, car 80 % des mails ont moins de 25 mots mean(sapply(1:length(x_token),function(x) length(x_token[[x]]))) # nombre de mots moyen par message quantile((sapply(1:length(x_token),function(x) length(x_token[[x]]))),0.8) # quantiles du nb de mots par message maxlength <- 25 x_train <- pad_sequences(x_encode, maxlen = maxlength) #x_train <- pad_sequences(x_encode, maxlen = maxlength, padding = "post", truncating = "post") dim(x_train) # effet des transformations successives head(x_token,1) head(x_encode,6) x_train[1,] # échantillon de test x_token <- tokenize_words(x_test, lowercase = FALSE, stopwords = NULL, strip_punct = TRUE, strip_numeric = FALSE) x_encode <- lapply(x_token, function(x) text_one_hot(paste(x, collapse=" "), maxwords)) x_test <- pad_sequences(x_encode, maxlen = maxlength) #x_test <- pad_sequences(x_encode, maxlen = maxlength, padding = "post", truncating = "post") dim(x_test) y_train <- to_categorical(y_train-1, num_classes) # on soustrait 1 car la variable initiale est 1/2 y_test <- to_categorical(y_test-1, num_classes) # définition du modèle LSTM model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 50, input_length=maxlength) %>% layer_lstm(units = 128, recurrent_dropout = 0.2) %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle GRU model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 50, input_length=maxlength) %>% layer_gru(units = 128, recurrent_dropout = 0.2) %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle GRU bidirectionnel model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 50, input_length=maxlength) %>% bidirectional(layer_gru(units = 128)) %>% #layer_dropout(0.2, seed=235) %>% layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # définition du modèle CNN model <- keras_model_sequential() model %>% layer_embedding(input_dim = maxwords, output_dim = 50, input_length=maxlength) %>% layer_conv_1d(filters = 128, kernel_size = 5) %>% layer_batch_normalization() %>% layer_activation("relu") %>% layer_global_max_pooling_1d() %>% # max pooling global ou aplatissement layer_dense(units = 64, activation = "relu") %>% layer_dropout(0.5, seed=235) %>% layer_dense(units = num_classes, activation = 'softmax') summary(model) # compilation model %>% compile(loss = 'categorical_crossentropy', optimizer = 'adam', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = 'adadelta', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = 'adagrad', metrics = c('accuracy')) model %>% compile(loss = 'categorical_crossentropy', optimizer = optimizer_rmsprop(lr=0.001), metrics = c('accuracy')) # apprentissage system.time(history <- fit(model, x_train, y_train, batch_size = 128, epochs = 10, validation_data = list(x_test, y_test))) # prédictions pred <- predict_classes(model, x_test) test.label <- max.col(y_test) - 1 table(test.label, pred) mean(test.label != pred) # historique de l'apprentissage plot(history) # Section 10.6 : Analyse de sentiments # ------------------------------------------------------------------- # version du CRAN install.packages("sentimenr") # version la plus récente remove.packages(pkgs = c("sentimentr")) library(devtools) install_github("trinker/sentimentr") library(sentimentr) # Incipit Pride and Prejudice txt <- c(" It is a truth universally acknowledged, that a single man in possession of a good fortune, must be in want of a wife. However little known the feelings or views of such a man may be on his first entering a neighbourhood, this truth is so well fixed in the minds of the surrounding families, that he is considered the rightful property of some one or other of their daughters. 'My dear Mr. Bennet,' said his lady to him one day, 'have you heard that Netherfield Park is let at last?' Mr. Bennet replied that he had not. 'But it is,' returned she; 'for Mrs. Long has just been here, and she told me all about it.' Mr. Bennet made no answer. 'Do you not want to know who has taken it?' cried his wife impatiently. 'You want to tell me, and I have no objection to hearing it.'") (sentim <- sentiment(txt)) (sentim <- sentiment_by(txt)) # coloriage des sentiments dans un fichier HTML highlight(sentim) # il faut avoir la version 2.8.0 du package sentimentr # phrases de test (code absent du livre) txt <- c("I like it", "I do not like it", "I don’t like it", "I really like it", "I dislike it", "I really dislike it") # exemples de textes data(hotel_reviews) txt <- hotel_reviews$text # sentiments des phrases sentim <- sentiment(txt) sentim$no <- seq(1,nrow(sentim)) head(sentim,30) # extraction des mots positifs et négatifs de chaque phrase pol_words <- extract_sentiment_terms(txt) head(pol_words) # phrases analysées head(pol_words$sentence) # phrase la plus positive pol_words$sentence[sentim[which.max(sentim$sentiment)]$no] # phrase la plus négative pol_words$sentence[sentim[which.min(sentim$sentiment)]$no] # évolution des sentiments au fil des textes sentimg <- sentiment_by(txt) head(sentimg) plot(sentimg$ave_sentiment, type="l", main="Exemple de trajectoire", xlab = "temps de la narration", ylab= "polarité") # classement des mots selon leur polarité mots <- attributes(pol_words)$counts mots <- mots[order(mots$polarity, decreasing=T), ] # mots les plus positifs head(mots,5) # mots les plus négatifs tail(mots,5) # polarité des mots de chaque phrase head(pol_words$neutral) head(pol_words$positive) tail(pol_words$negative) # --------------------------------------------------------------------------------------------------------- # Chapitre 11 : Analyse des réseaux sociaux # --------------------------------------------------------------------------------------------------------- # Section 11.5 : Graphes avec R # ------------------------------------------------------------------- library(igraph) # exemples simples graphe <- graph.formula( Alice-Bob-Cecil-Alice, Daniel-Cecil-Eugene, Cecil-Gordon ) graphe <- graph.formula( Alice +-+ Bob --+ Cecil +-- Daniel, Eugene --+ Gordon:Helen ) plot(graphe) # coloriage des arêtes E(graphe)$color <- "yellow" # taille et couleur des sommets V(graphe)$size <- 20 V(graphe)$color <- "red" plot(graphe) plot(graphe, vertex.label=V(graphe)$name, vertex.size=1, vertex.label.cex=2, edge.arrow.size=1, layout= layout.fruchterman.reingold) tkplot(graphe) graphe <- graph_from_data_frame(liens, vertices=acteurs, directed=TRUE) graphe <- graph_from_data_frame(liens, vertices=NULL, directed=TRUE) # analyse du graphe summary(graphe) V(graphe) # liste des sommets E(graphe) # liste des arêtes V(graphe)$name # noms des sommets length(V(graphe)) # nombre de sommets vcount(graphe) length(E(graphe)) # nombre d'arêtes ecount(graphe) edge_attr(graphe) # attributs du graphe table(E(graphe)$type) # composantes connexes is.connected(graphe) # le graphe est-il connexe no.clusters(graphe) # nb de composantes connexes clusters(graphe) # taille des composantes et appartenance sommets table(clusters(graphe)$csize) # distribution de la taille des composantes # réduction du graphe à sa première composante connexe vertices <- which(clusters(graphe)$membership==1) graphe <- induced.subgraph(graphe,vertices) # degré d'un noeud = nb d'arêtes adjacentes table(degree(graphe)) # distribution du degré des sommets table(degree(graphe,mode="in")) # degré des sommets entrants table(degree(graphe,mode="out")) # degré des sommets sortants max(degree(graphe)) max(degree(graphe, mode="in")) max(degree(graphe, mode="out")) plot(degree.distribution(graphe, mode="in"), log="xy") plot(degree.distribution(graphe), log="xy") which.max(degree(graphe)) # sommet ayant le plus d'arêtes which.max(degree(graphe, mode="in")) # sommet ayant le plus d'arêtes entrantes which.max(degree(graphe, mode="out")) # sommet ayant le plus d'arêtes sortantes which(degree(graphe, mode="in")==max(degree(graphe, mode="in"))) # sommets ayant le plus d'arêtes entrantes which(degree(graphe, mode="out")==max(degree(graphe, mode="out"))) # sommets ayant le plus d'arêtes sortantes which(degree(graphe)==1) # sommets ayant un seul voisin # mesures d'influence graph.density(graphe) # densité du graphe transitivity(graphe) # coefficient de clustering summary(betweenness(graphe)) # intermédiarité plot(closeness(graphe)) # centralité de proximité summary(eccentricity(graphe)) # excentricité des sommets # diamètre diameter(graphe, directed=F) diameter(graphe, directed=T) # distance en prenant les chemins avec les arcs à l’endroit seulement radius(graphe) # Section 11.5 : Graphes de packages # ------------------------------------------------------------------- library(miniCRAN) pkgDep("miniCRAN", suggests = FALSE, enhances = FALSE) # une sélection de 15 packages du CRAN packages <- c("igraph", "ggplot2", "data.table", "glmnet", "shiny", "nnet", "randomForest", "rpart", "xgboost", "mxnet", "keras", "ff", "ffbase", "h2o", "tm") graphe <- makeDepGraph(packages, suggests = FALSE, enhances = FALSE) plot(graphe, cex=1.25, legendPosition = c(1, -1)) library(igraph) # stats de base sur le graphe length(V(graphe)) # nombre de sommets length(E(graphe)) # nombre d'arêtes table(E(graphe)$type) table(degree(graphe)) # distribution du degré des sommets table(degree(graphe, mode="in")) # degré des sommets entrants table(degree(graphe, mode="out")) which(degree(graphe, mode="in")==max(degree(graphe, mode="in"))) # sommets ayant le plus d'arêtes entrantes which(degree(graphe, mode="out")==max(degree(graphe, mode="out"))) # sommets ayant le plus d'arêtes sortantes # graphe sans échelle ? fit <- fit_power_law(degree(graphe)) fit # composantes connexes is.connected(graphe) no.clusters(graphe) table(clusters(graphe)$csize) table(clusters(graphe)$csize)/length(V(graphe)) # indicateurs graph.density(graphe) transitivity(graphe) # diamètre diameter(graphe, directed=F) diameter(graphe, directed=T) # graphe petit-monde ? mean_distance(graphe) log(length(V(graphe))) # pageRank page_rank(graphe)$vector summary(page_rank(graphe)$vector) quantile(page_rank(graphe)$vector,0.9) names(which(page_rank(graphe)$vector > quantile(page_rank(graphe)$vector,0.9))) which(page_rank(graphe)$vector==max(page_rank(graphe)$vector)) #save(graphe, file="D:/Data Mining/Livres Data Mining/Technip Big Data/graphe_15packages.RData") #load("D:/Data Mining/Livres Data Mining/Technip Big Data/graphe_15packages.RData") # coloriage des sommets à fort pagerank V(graphe)$size <- 2 V(graphe)$size[which(page_rank(graphe)$vector > quantile(page_rank(graphe)$vector,0.9))] <- 10 V(graphe)$color <- "white" V(graphe)$color[which(page_rank(graphe)$vector > quantile(page_rank(graphe)$vector,0.9))] <- "red" # coloriage des arêtes E(graphe)$color <- "grey" col <- rep("blue", length(E(graphe)$type)) col[E(graphe)$type=="Imports"] <- "red" col[E(graphe)$type=="Depends"] <- "yellow" col[E(graphe)$type=="LinkingTo"] <- "black" plot.igraph(graphe, layout=layout.fruchterman.reingold) tkplot(graphe, layout=layout.fruchterman.reingold, vertex.color=V(graphe)$color, vertex.size=V(graphe)$size, edge.color=col, vertex.label.cex=2) # tous les packages du CRAN packages <- row.names(available.packages()) graphe <- makeDepGraph(packages, suggests = FALSE, enhances = FALSE) max(degree(graphe, mode="in")) which(degree(graphe, mode="in")==max(degree(graphe, mode="in"))) max(degree(graphe, mode="out")) which(degree(graphe, mode="out")==max(degree(graphe, mode="out"))) #save(graphe, file="D:/Data Mining/Livres Data Mining/Technip Big Data/graphe_tous_packages.RData") #load("D:/Data Mining/Livres Data Mining/Technip Big Data/graphe_tous_packages.RData") # graphe des frontières de pays liens <- read.csv("Frontieres_pays.csv") liens <- read.csv2("D:/Data Mining/Logiciel R/Programmes R/Graphes/Frontieres_pays.csv") graphe <- graph.data.frame(liens, vertices=NULL, directed=TRUE) summary(graphe) graphe <- as.undirected(graphe) vcount(graphe) ecount(graphe) table(degree(graphe)) which.max(degree(graphe)) # Section 11.6.5 : Détection des communautés avec R # ------------------------------------------------------------------- # détection des communautés cm <- fastgreedy.community(as.undirected(graphe)) # KO si arêtes multiples cm <- fastgreedy.community(simplify(graphe)) # KO si graphe orienté system.time(cm <- fastgreedy.community(simplify(as.undirected(graphe)))) # OK system.time(cm <- edge.betweenness.community(graphe)) # KO si graphe non orienté system.time(cm <- edge.betweenness.community(as.undirected(graphe))) # calcul très long sur un grand graphe system.time(cm <- cluster_walktrap(graphe)) system.time(cm <- walktrap.community(graphe)) system.time(cm <- cluster_louvain(graphe)) # KO si graphe orienté system.time(cm <- cluster_louvain(as.undirected(graphe))) system.time(cm <- cluster_leading_eigen(graphe)) # warning mais pas KO si graphe orienté system.time(cm <- cluster_leading_eigen(as.undirected(graphe))) system.time(cm <- cluster_infomap(graphe)) # warning mais pas KO si graphe orienté system.time(cm <- cluster_infomap(as.undirected(graphe))) system.time(cm <- cluster_label_prop(graphe)) # calcul très long sur un grand graphe system.time(cm <- cluster_spinglass(graphe, spins=4)) # KO si graphe non connexe system.time(cm <- cluster_spinglass(graphe, vertex=1)) # KO si graphe non connexe # réduction du graphe à sa première composante connexe vertices <- which(clusters(graphe)$membership==1) graphe <- induced.subgraph(graphe,vertices) # autre package : modMax library(modMax) # matrice d'adjacence adj <- get.adjacency(as.undirected(graphe)) # le graphe ne doit pas être orienté system.time(cm <- simulatedAnnealing(adj, fixed=10)) # calcul long sur un grand graphe system.time(cm <- louvain(adj)) # beaucoup plus long qu'avec le package igraph # méthode greedy standard system.time(cm <- greedy(adj)) # besoin de beaucoup de RAM et très long => utiliser fastgreedy du package igraph system.time(cm <- greedy(adj, refine = "fast")) # encore beaucoup plus long qu'avec le package igraph # Wakita, K., and T. Tsurumi system.time(cm <- greedy(adj, q = "wakita3")) # besoin de beaucoup de RAM pour wakita1, wakita2 et wakita3 commun <- cbind(cm[[3]], V(graphe)) commun[order(commun[,1]),] # sommets de chaque communauté # autre package : MCL library(MCL) # Markov Cluster Algorithm (MCA) de van Dongen adj <- get.adjacency(graphe) # le graphe peut être orienté system.time(cm <- mcl(adj, addLoops = TRUE, expansion = 2, inflation = 2, allow1 = FALSE, max.iter = 100)) # autre méthode spectrale (code absent du livre) #install.packages('motifcluster') library(motifcluster) #motif_cluster <- run_motif_clustering(adj, motif_name = "M1", motif_type = "func", mam_weight_type = "mean", mam_method = "sparse", type_lap = "rw", num_eigs = 4, num_clusts = 3) adj <- get.adjacency(graphe) # KO adj <- get.adjacency(as.undirected(graphe)) # OK system.time(motif_cluster <- run_motif_clustering(adj, motif_name = "M1", motif_type = "func", type_lap = "rw", num_eigs = 8, num_clusts = 8, restrict = TRUE)) motif_cluster$clusts length(motif_cluster$clusts) table(motif_cluster$clusts) motif_cluster$comps commun <- cbind(motif_cluster$clusts, V(graphe)[motif_cluster$comps]) commun[order(commun[,1]),] # sommets de chaque communauté names(V(graphe)) row.names(commun) commun cm <- merge(data.frame(V=names(V(graphe))), data.frame(V=row.names(commun), seg=commun[,1]), all.x=TRUE, sort = TRUE) cm$seg[is.na(cm$seg)] <- 0 cm # caractéristiques des communautés cm length(cm) # nombre de communautés sizes(cm) # taille des communautés table(sizes(cm)) # distribution de la taille des communautés head(membership(cm),20) # coloriage des communautés (package igraph) colbar <- rainbow(length(cm)) col <- colbar[membership(cm)] # coloriage des communautés (packages modMax et MCL) colbar <- rainbow(cm[[1]]) col <- colbar[cm[[3]]] # graphe statique layout1 <- layout.fruchterman.reingold plot(graphe, vertex.color=col, vertex.size=8, vertex.label.cex=0.85, layout=layout1) # graphe interactif tkplot(graphe, vertex.color=col) # Sections 11.10.1 et 11.10.4 : Collecte des tweets # ------------------------------------------------------------------- # création d'une application Twitter Aller à https://dev.twitter.com/ # les infos suivantes sont récupérées sur le compte Twitter de l'application : consumer_key <- "YOUR CONSUMER KEY" consumer_secret <- "YOUR CONSUMER SECRET" access_token <- "YOUR ACCESS TOKEN" access_secret <- "YOUR ACCESS SECRET" # authentification sur Twitter library(twitteR) setup_twitter_oauth(consumer_key, consumer_secret, access_token, access_secret) # autre package library(rtweet) app <- "nom_application" twitter_token <- create_token(app, consumer_key, consumer_secret, access_token, access_secret) # requêtes sur utilisateurs util <- getUser('@Pontifex_fr') util$followersCount # nombre de followers util$friendsCount # nombre de followees (abonnements) util$getFriends(n=10) util$getFavorites(n=10) util$getFollowers(n=10) userTimeline('@Pontifex_fr', n=100) # fil d’un utilisateur # tendances trend <- availableTrendLocations() # tendances head(trend) getTrends(trend[which(trend$country=="France"),]$woeid) # tendances sur la France getTrends(1) # tendances sur le Monde # extraction de tweets avec le package twitteR tweets <- searchTwitter("#COVID", lang='fr', n=15000, since='2020-06-01', until='2020-06-30') # extraction de tweets avec le package rtweet tweets <- search_tweets("#IA", "lang:fr", n = 10000, include_rts = FALSE, retryonratelimit = TRUE) # transformation de la liste des tweets en data frame system.time(df <- do.call("rbind", lapply(tweets, as.data.frame))) # commande directe et équivalente system.time(df <- twListToDF(tweets)) # autre commande un peu plus rapide library(data.table) system.time(df <- rbindlist(lapply(tweets, as.data.frame))) head(df) # dates des tweets table(as.Date(df$created)) # auteurs des tweets length(unique(df$screenName)) # nombre d'auteurs # auteurs les plus fréquents d <- as.data.frame(table(df$screenName)) d <- d[order(d$Freq, decreasing=T), ] names(d) <- c("User","Tweets") head(d,10) # autres informations table(df$longitude) # localisation table(df$latitude) table(df$isRetweet) # messages qui sont des re-tweets (T/F) length(which(!df$isRetweet)) # nombre de tweets originaux length(which(df$isRetweet)) # nombre de retweets table(df$retweeted) # messages retweetés (T/F) table(df$retweetCount) # nombre de re-tweets du tweet original summary(df$retweetCount) # message le plus retweeté df[df$retweetCount==max(df$retweetCount), c('created','screenName','text','retweetCount')] # auteurs les plus fréquents des tweets qui ne sont pas des retweets d <- as.data.frame(table(df$screenName[which(!df$isRetweet)])) d <- d[order(d$Freq, decreasing=T), ] names(d) <- c("User","Tweets") head(d,10) # suppression des doublons (à cause des retweets multiples) dim(df) df <- df[!duplicated(df$text),] dim(df) # source des tweets d <- as.data.frame(table(df$statusSource)) d <- d[order(d$Freq, decreasing=T), ] names(d) <- c("Source", "Tweets") head(d) # analyse des hashtags # récupérer l'ensemble des mots (délimités par des ESPACE) mots <- unlist(strsplit(df$text," ")) # un mot est un hashtag s'il débute par # is_hashtag <- regexpr("^#[[:alnum:]_]*", mots) # récupérer l'ensemble des thèmes désignés par un "#" dans les messages liste_hashtags <- regmatches(mots, is_hashtag) # nombre de hashtags recensés length(liste_hashtags) length(unique(liste_hashtags)) # nombre d'apparition de chaque hashtag trié selon la fréquence décroissante hashtags <- sort(table(liste_hashtags), decreasing=TRUE) # affichage des 10 hashtags les plus fréquents head(hashtags, 10) #liste des hashtags contenant IA et leur nombre d'appartition nb_hashtags <- table(liste_hashtags) hashtags_ia <- nb_hashtags[grep("IA", names(nb_hashtags))] print(sort(hashtags_ia, decreasing=TRUE)) # wordcloud des hashtags library(wordcloud) wordcloud(names(hashtags)[-1], hashtags[-1], scale=c(3,.5), random.order = FALSE, colors=brewer.pal(6, "Dark2")) # exclusion des trois hashtags les plus fréquents wordcloud(names(hashtags)[-(1:3)], hashtags[-(1:3)], scale=c(3,.5), min.freq=10, random.order = FALSE, rot.per=.2, max.words=1000) # Sections 11.10.2 et 11.10.4 (2e partie) : Mise en forme du corpus # ------------------------------------------------------------------- # suppression des caractères non imprimables texte <- unlist(lapply(df$text, function(x) gsub("[^[:print:]]", "", x))) # création du corpus library(tm) myCorpus <- Corpus(VectorSource(texte)) # affichage du contenu des documents inspect(myCorpus[1:6]) length(unlist(myCorpus)) # conversion de tous les caractères en bas de casse myCorpus <- tm_map(myCorpus, content_transformer(tolower)) #myCorpus <- tm_map(myCorpus, tolower) # suppression des mots commençant par @ (avant suppression ponctuation !) myCorpus <- tm_map(myCorpus, function(x) gsub("@[[:alnum:]]*", "", x)) #myCorpus <- tm_map(myCorpus, function(x) gsub("@\\w+ *", "", x)) # variante # suppression des hashtags # (à faire au début des transformations) myCorpus <- tm_map(myCorpus, function(x) gsub("#", " ", x)) # suppression des URL (avant avoir enlevé la ponctuation :) myCorpus <- tm_map(myCorpus, function(x) gsub("http[[:graph:]]*", "", x)) # suppression des nombres #myCorpus <- tm_map(myCorpus, removeNumbers) # définition des mots vides (stopwords), en y ajoutant des mots personnels stopwords('french') myStopwords <- c(stopwords('french'), "via","amp","a","rt","où","ça") myCorpus <- tm_map(myCorpus, removeWords, myStopwords) # suppression de la ponctuation (sauf dans les mots composés) myCorpus <- tm_map(myCorpus, function(x) gsub("[[:punct:]]* *(\\w+[&'-]\\w+)|[[:punct:]]+ *| {2,}", " \\1", x)) # suppression de la ponctuation (sauf les - dans les mots composés) #myCorpus <- tm_map(myCorpus, function(x) removePunctuation(x, preserve_intra_word_dashes = TRUE)) # suppression des blancs inutiles myCorpus <- tm_map(myCorpus, stripWhitespace) # passage d'un corpus à un texte texte.corpus <- data.frame(text=head(unlist(myCorpus),length(texte)), date=as.Date(df$created), lon=df$longitude, lat=df$latitude, stringsAsFactors=F) head(texte.corpus) # Section 11.10.4 (2e partie) : Exemple # ------------------------------------------------------------------- # remplacement des abréviations myCorpus <- tm_map(myCorpus, function(x) gsub("\\", "intelligence artificielle", x)) # correction d'erreurs myCorpus <- tm_map(myCorpus, gsub, pattern="covid-19", replacement="covid19") myCorpus <- tm_map(myCorpus, gsub, pattern="covid_19", replacement="covid19") myCorpus <- tm_map(myCorpus, gsub, pattern="lintelligence", replacement="intelligence") myCorpus <- tm_map(myCorpus, gsub, pattern="dintelligence", replacement="intelligence") myCorpus <- tm_map(myCorpus, gsub, pattern="intelligenceartificielle", replacement="intelligence artificielle") myCorpus <- tm_map(myCorpus, gsub, pattern="machinelearning", replacement="machine learning") myCorpus <- tm_map(myCorpus, gsub, pattern="deeplearning", replacement="deep learning") myCorpus <- tm_map(myCorpus, gsub, pattern="marketingdigital", replacement="marketing digital") # suppression des mots courts myCorpus <- tm_map(myCorpus, function(x) gsub("\\<\\w{1,2}\\>", "", x)) # lemmatisation en ne conservant que les mots trouvés dans le dictionnaire et qui sont des NOMS ou des ADJECTIFS ou des VERBES # NB : à lancer avant la suppression des accents, des ' et le regroupement des expressions, mais après la suppression des @ et # et corrections des erreurs library(data.table) lemme <- fread("D:/Data Mining/Logiciel R/Programmes R/Text Mining/Lemmatisation/lemmes-fr.txt", sep="\t", header=F) setkey(lemme,V1) head(lemme);tail(lemme) myCorpus <- tm_map(myCorpus, function(y) unlist(lapply(lapply(strsplit(y,' '), function(x) as.character(ifelse(!is.na(lemme[x]$V2)&(lemme[x]$V3 %in% c("NOM", "ADJ")), lemme[x]$V2, ""))),function(x) paste(x,collapse =" ")))) inspect(myCorpus[1:10]) # suppression des accents myCorpus <- tm_map(myCorpus, function(x) gsub("[éèêÉÈÊ]","e", x)) myCorpus <- tm_map(myCorpus, function(x) gsub("[àâÀÂÂ]","a", x)) myCorpus <- tm_map(myCorpus, function(x) gsub("[ùûüÙÛÜ]","u", x)) myCorpus <- tm_map(myCorpus, function(x) gsub("[îïÎÏ]","i", x)) myCorpus <- tm_map(myCorpus, function(x) gsub("[ôÔ]","o", x)) myCorpus <- tm_map(myCorpus, function(x) gsub("ç","c", x)) # suppression des blancs inutiles myCorpus <- tm_map(myCorpus, stripWhitespace) # affichage du corpus inspect(myCorpus[1:10]) # matrice tdm tdm <- TermDocumentMatrix(myCorpus, control = list(wordLengths=c(3,Inf))) as.matrix(tdm[1:10,1:10]) tdm # fréquence des termes Zipf_plot(tdm) findFreqTerms(tdm, lowfreq=100) freq <- sort(rowSums(as.matrix(tdm)), decreasing=TRUE) d <- data.frame(freq) head(d,15) library(wordcloud) wordcloud(rownames(d), d$freq, scale=c(3,.7), min.freq=5, random.order = FALSE, rot.per=.2, max.words=100) # Section 11.10.5 : Classification de termes et de documents # ------------------------------------------------------------------- # CAH tdm2 <- removeSparseTerms(tdm, sparse=0.985) m2 <- as.matrix(tdm2) distMatrix <- dist(scale(m2)) h <- hclust(distMatrix, method="ward.D2") plot(h) # calcul inertie interclasse et totale : à la main # inertie totale totss <- function(dmatrix) { globalmean <- apply(dmatrix, 2, mean) sum(apply(dmatrix, 1, function(x) {(sum((x-globalmean)^2, na.rm = T))})) } (totalss <- totss(scale(m2))) # inertie intra-classe totale wss.total <- function(dmatrix, groupes) { wsstot <- 0 l <- length(unique(groupes)) for(i in 1:l) wsstot <- wsstot + totss(subset(dmatrix, groupes==i)) wsstot } # valeurs d'inertie intra-classe entre 1 et kmax classes kmax <- 20 wss <- numeric(kmax) for(k in 1:kmax) { labels <- cutree(h, k=k) wss[k] <- wss.total(scale(m2), labels) } # inertie interclasse bss <- totalss - wss # calcul du R² (ev = explained variance) et de l'indice de Calinski-Harabasz pour chaque nb de classes possible n <- dim(scale(m2))[1] (inertie <- cbind(k=1:kmax, wss, bss, ev=bss/totalss, ch=(bss/0:(kmax-1))/(wss/n-1:kmax))) # maximisation du saut de R² semi partiel (spr2 <- inertie[2:20,4 ] - inertie[1:19, 4]) plot(spr2, type="b", pch=16) sapply(1:19, function(n) (spr2[n]/spr2[n+1])) 1+which.max(sapply(1:19, function(n) (spr2[n]/spr2[n+1]))) rank(sapply(1:19, function(n) (spr2[n]/spr2[n+1])), na.last = NA) # indice de Calinski-Harabasz tdm2 <- removeSparseTerms(tdm, sparse=0.99) m2 <- as.matrix(tdm2) #m2 <- m2[, sample(ncol(m2), 15000)] dim(m2) m3 <- t(m2) library(fpc) clustering.ch <- kmeansruns(m3, krange=2:8, criterion="ch") clustering.ch$crit (k <- clustering.ch$bestk) # k-means set.seed(123) kmeansResult <- kmeans(m3, k, nstart=10) # top des termes dans chaque classe for (i in 1:k) { cat(paste("cluster ", i, " : ", sep="")) s <- sort(kmeansResult$centers[i,], decreasing=T) cat(names(s)[1:5], "\n") } # taille des classes kmeansResult$size # Section 11.10.6 : Score d’opinion # ------------------------------------------------------------------- # fonction renvoyant le "score de sentiment" d'une phrase library(stringr) f <- function(sentence) { # récupération des mots de la phrase dans un vecteur # str_split est plus rapide que strsplit words = unlist(str_split(sentence, '\\s+')) # recherche des mots du vecteur dans la liste des mots positifs / négatifs # la fonction "match" renvoie la position du mots du vecteur dans la liste des mots # mais c'est juste sa présence dans la liste qui nous intéresse pos.matches = !is.na(match(words, pos.words)) neg.matches = !is.na(match(words, neg.words)) # le score est le nombre de mots (dans la liste des) positifs - le nombre de mots négatifs score = sum(pos.matches) - sum(neg.matches) return(score) } # utilisation de la base FEEL words <- read.csv2("D:/Data/NLP/FEEL.csv") pos.words <- as.vector(words[words$polarity=="positive","word"]) neg.words <- as.vector(words[words$polarity=="negative","word"]) Encoding(pos.words) <- "UTF-8" # pour conserver les accents sans les remplacer par des caractères spéciaux Encoding(neg.words) <- "UTF-8" # pour conserver les accents sans les remplacer par des caractères spéciaux # calcul du score des textes du corpus texte.corpus <- data.frame(text=unlist(myCorpus)[1:length(myCorpus)], date=as.Date(df$created), stringsAsFactors=F) head(texte.corpus) score <- Vectorize(f)(unlist(texte.corpus$text)) table(score) # évolution du score d'opinion dans le temps evol.score <- aggregate(score, by=list(texte.corpus$date), summary) plot(evol.score$Group.1, evol.score$x[,4], type="l", xlab="Date", ylab="", ylim=c(floor(range(evol.score$x[,4])[1]),ceiling(range(evol.score$x[,4])[2])), col="red", lab=c(10,10,0), cex.axis=0.1) axis(2, cex.axis=1.5) lines(evol.score$Group.1, evol.score$x[,3], type="l", xlab="Date", ylab="", lty=3, col="blue") mtext("Score moyen", 2, line=2, col="red", cex=2) mtext("Score médian", 4, line=0.5, col="blue", cex=2) axis.Date(side=1, cex.axis=1.2, at=seq(range(evol.score$Group.1)[1], range(evol.score$Group.1)[2], "days")) # Section 11.10.7 : Graphe des termes avec leur connotation # ------------------------------------------------------------------- # pas plus de 99% de termes nuls dans la matrice DTM tdm2 <- removeSparseTerms(tdm, sparse=0.99) m2 <- as.matrix(tdm2) m2[1:10,1:20] # transformation en matrice booléenne m2[m2>=1] <- 1 # transformation en matrice d'adjacence termMatrix <- m2 %*% t(m2) termMatrix[1:20,1:20] # graphe de co-occurrences des termes dans des documents library(igraph) # deux termes sont reliés s'ils sont co-occurents dans un document g <- graph.adjacency(termMatrix, weighted=T, mode="undirected") # suppression des boucles et arêtes multiples g <- simplify(g) # set labels and degrees of vertices V(g)$label <- V(g)$name V(g)$degree <- degree(g) # taille des sommets en fonction de la fréquence des termes freq <- sort(rowSums(m2), decreasing=TRUE) d <- data.frame(word=names(freq), freq=freq) V(g)$size <- log(d[match(V(g)$name,d$word),2])*5 V(g)$color <- "grey" # affichage du graphe d’adjacence des termes set.seed(123) plot(g, vertex.label=V(g)$name, vertex.size=V(g)$size, vertex.label.cex=1, edge.arrow.size=0.5, layout=layout.fruchterman.reingold) # graphe interactif library(tcltk) tkplot(g) # détection des communautés cm <- fastgreedy.community(g) cm <- edge.betweenness.community(g) cm <- walktrap.community(g) cm <- cluster_louvain(g) length(cm) # nombre de communautés sizes(cm) # taille des communautés # coloriage des communautés colbar <- rainbow(length(cm)) col <- colbar[membership(cm)] # coloriage en gris library(RColorBrewer) colbar <- brewer.pal(length(cm), "Greys") col <- colbar[membership(cm)] # affichage des communautés du graphe d’adjacence set.seed(123) plot(g, vertex.color=col, layout=layout.fruchterman.reingold) plot(g, mark.groups=communities(cm), vertex.color=col, layout=layout.fruchterman.reingold) # graphe interactif tkplot(g, vertex.color=col)