diff --git a/code/ODL.R b/code/ODL.R
new file mode 100755
index 0000000000000000000000000000000000000000..99a3c31694bec2762522fd44c3ad9b0eb90f8460
--- /dev/null
+++ b/code/ODL.R
@@ -0,0 +1,502 @@
+rm(list=ls())
+require(assertthat)
+require(lars)
+require(ggplot2)
+require(corrplot)
+require(ggcorrplot)
+#setwd("/Volumes/TCDATA/Projets/MNIST")
+setwd("~/Documents/Cours/PRR/MNIST")
+source("./skeleton.R")
+
+# ----------------------------------------------------------
+# Plot permettant de visualiser les features de X (mnist$test$x)
+plot_mnist_random()
+
+# ----------------------------------------------------------
+# norm_dico
+#
+# Require :
+# D := un dictionnaire
+#
+# Ensure : Retourne D normalisé pour que D appartienne à C
+# ----------------------------------------------------------
+norm_dico <- function(D) {
+ dim_D <- dim(D) ; p = dim_D[1] ; k = dim_D[2]
+ for(j in 1:k) {
+ dj <- D[,j]
+ D[,j] <- dj / max(1,sqrt(sum(dj^2)))
+ }
+ D
+}
+
+# ----------------------------------------------------------
+# make_dico_1
+#
+# Require :
+# p := nombre de features du dictionnaire
+# k := nombre d'atoms du dictionnaire
+#
+# Ensure : Retourne un dictionnaire de k atoms exprimés
+# par p parameters comme les features de X dont nous
+# souhaitons construire un dictionnaire le représentant
+#
+# Ensure+ : Les atoms sont aléatoires
+# ----------------------------------------------------------
+make_dico_1 <- function(p,k) {
+ norm_dico(matrix(data = rnorm(p*k), nrow = p, ncol = k))
+}
+
+plot_mnist_choose(t(make_dico_1(784,100)),10)
+
+# ----------------------------------------------------------
+# make_dico_2
+#
+# Require :
+# k := nombre d'atoms du dictionnaire
+#
+# Ensure : Retourne un dictionnaire de k atoms exprimés
+# par p parameters comme les features de X dont nous
+# souhaitons construire un dictionnaire le représentant
+#
+# Ensure+ : Les k atoms sont pris dans X = mnist$train$x
+# et p = 784
+# ----------------------------------------------------------
+make_dico_2 <- function(k) {
+ dim_X <- dim(mnist$train$x) ; n = dim_X[1] ; p = dim_X[2]
+ norm_dico(t(mnist$train$x[sample(1:n,k),]))
+}
+
+plot_mnist_choose(t(make_dico_2(100)),10)
+
+# ----------------------------------------------------------
+# make_dico_3
+#
+# Require :
+# m := nombre de features du dictionnaire
+# k := nombre d'atoms du dictionnaire
+#
+# Ensure : Retourne un dictionnaire de k atoms exprimés
+# par m parameters comme les features de X dont nous
+# souhaitons construire un dictionnaire le représentant
+#
+# Ensure+ : Les k atoms ont un schéma carré de 9px d'aire
+# situé aléatoirement
+# ----------------------------------------------------------
+make_dico_3 <- function(m,k) {
+ D <- matrix(data = 0, nrow = m, ncol = k)
+ p <- sqrt(m)
+ for(t in 1:k) {
+ i = floor(runif(1,2,p))-1
+ j = floor(runif(1,2,p))
+ D[1+i-1+p*(j-1),t] <- 1 ; D[1+i+p*(j-1),t] <- 1 ; D[1+i+1+p*(j-1),t] <- 1
+ D[1+i-1+p*(j),t] <- 1 ; D[1+i+p*(j),t] <- 1 ; D[1+i+1+p*(j),t] <- 1
+ D[1+i-1+p*(j-2),t] <- 1 ; D[1+i+p*(j-2),t] <- 1 ; D[1+i+1+p*(j-2),t] <- 1
+ }
+ norm_dico(D)
+}
+
+plot_mnist_choose(t(make_dico_3(784,100)),10)
+
+# ----------------------------------------------------------
+# make_dico_4
+#
+# Require :
+# k := 10k est le nombre d'atoms du dictionnaire
+#
+# Ensure : Retourne un dictionnaire de k atoms exprimés
+# par p parameters comme les features de X dont nous
+# souhaitons construire un dictionnaire le représentant
+#
+# Ensure+ : Les 10k atoms sont pris dans X = mnist$train$x
+# et p = 784
+# Il y a k atoms de chaque chiffre de 0 à 9
+# ----------------------------------------------------------
+make_dico_4 <- function(k) {
+ D <- c()
+ for(i in 0:9) {
+ ind <- sample(which(mnist$test$y==i),k)
+ D <- cbind(D,t(mnist$test$x[ind,]))
+ }
+ norm_dico(D)
+}
+
+plot_mnist_choose(t(make_dico_4(10)),10)
+
+make_test <- function() { # like make_dico_4(1)
+ D <- c()
+ for(i in 0:9) {
+ ind <- sample(which(mnist$test$y==i),1)
+ D <- cbind(D,matrix(data=mnist$test$x[ind,], nrow = 784, ncol = 1))
+ }
+ norm_dico(D)
+}
+
+plot_mnist_choose(t(make_test()),5)
+
+# ----------------------------------------------------------
+# cost_function
+#
+# Require :
+#
+# Ensure :
+# ----------------------------------------------------------
+cost_function <- function(X,D,l) {
+
+ fn <- 0
+ dim_X <- dim(X) ; n <- dim_X[1] ; N <- n
+ norm_vec <- function(x) sqrt(sum(x^2))
+ cat("0 0.1 0.2 0.3 0.4 ")
+ cat("0.5 0.6 0.7 0.8 0.9 1\n")
+ cat ("|")
+
+ for(t in 1:n) {
+ ################################################################################
+ # Extraction d'une feature xt de X
+ xt <- X[t,]
+
+ ################################################################################
+ # Calcul de la décomposition rt de xt
+ try(model <- lars(D, xt, type = "lasso"), silent = T)
+ if(is.na(model)) {
+ N <- N-1
+ next
+ }
+ ind <- which.min(model$lambda>=l) # max(which(model$lambda>=l))
+ if(is.na(ind)) ind <- model$lambda[1]
+ rt <- coef(model)[ind,]
+
+ ################################################################################
+ # Calcul de la représentation yt de xt dans D
+ yt <- D%*%rt
+
+ ################################################################################
+ # Contribution de xt dans fn
+ fn <- fn + 0.5*norm_vec(xt-yt)^2 + l*norm_vec(rt)
+
+ # Visualisation de l'avancement de l'apprentissage
+ if(n>=100 & t %% floor(n/100) == 0) cat (".") # cat("Etape t =",t,"\n")
+ if(n>=10 & t %% floor(n/10) == 0) cat ("|")
+ }
+ cat("\n")
+
+ fn/N
+}
+
+# ----------------------------------------------------------
+# ----------------------------------------------------------
+# ODL := Online Dictionary Learning Algorithm
+#
+# Require :
+# X := train dataset (matrix : size = n*p)
+# lambda := regularization parameter (double : > 0)
+# D := initial dictionary (matrix : size = p*k)
+# Tt := update iteration number (int : > 0)
+# Tb := batch size (int : > 0)
+# Tc := converge iteration number (int : > 0)
+# ----------------------------------------------------------
+# ----------------------------------------------------------
+
+ODL <- function(X, lambda, D, Tt, Tb, Tc, alea=T) {
+
+ ################################################################################
+ # Initialisation des paramètres n, p, k
+
+ # n := nombre de features de X
+ # p := nombre de parameters de X
+ dim_X <- dim(X) ; n <- dim_X[1] ; p <- dim_X[2]
+ # k := nombre d'atoms de D
+ # p := nombre de parameters des atoms de D
+ dim_D <- dim(D) ; p2 <- dim_D[1] ; k <- dim(D)[2]
+
+ assert_that(p==p2, msg = "Dimension of X and D not compatible")
+
+ mu <- 1
+
+ ################################################################################
+ # Initialisation des matrices d'informations pour le block-coordinate descent
+ # At := sum(rt*rt^T)
+ At <- matrix(data = 0, nrow = k, ncol = k)
+ # Bt := sum(xt*rt^T)
+ Bt <- matrix(data = 0, nrow = p, ncol = k)
+
+ ################################################################################
+ # Initialisation de nos variables
+ # Dt := dictionary après t updates
+ Dt <- D
+ # xt := t-th feature de X pour la t-th update
+ xt <- c()
+ # rt := sparse coding de la t-th feature (xt)
+ rt <- c()
+ # cc := compteur
+ cc <- 0
+ # ti := indice de la t-th feature de X
+ ti <- 0
+ # ind_ti := ensemble des indices ti
+ ind_ti <- c()
+
+ cat("0 0.1 0.2 0.3 0.4 ")
+ cat("0.5 0.6 0.7 0.8 0.9 1\n")
+
+ ################################################################################
+ # Boucle pour le Online Dictionary Learning - t-ème mise à jour
+ cat ("|")
+ for(t in 1:Tt) {
+ ################################################################################
+ # Extraction d'une feature xt de X
+ if(alea) ti <- floor(runif(1,1,n+1))
+ else ti <- t
+ ind_ti <- c(ind_ti,ti)
+ xt <- X[ti,]
+
+ ################################################################################
+ # Sparse Coding Problem : Résolution par la méthode du LARS-Lasso
+ try(model <- lars(Dt, xt, type = "lasso"), silent = T)
+ if(is.na(model)) {
+ cc <- cc - 1
+ t <- t - 1
+ next
+ }
+ ind <- which.min(model$lambda>=lambda) # max(which(model$lambda>=lambda))
+ #if(is.na(ind)) {ind <- model$lambda[length(SC$lambda)]}
+ assert_that(!is.na(ind), msg = "Trouble with value of lambda")
+ rt <- coef(model)[ind,]
+
+ ################################################################################
+ # Mise à jour des matrices d'informations
+ At <- mu*At + matrix(data = rt, ncol = 1) %*% t(matrix(data = rt, ncol = 1))
+ Bt <- mu*Bt + matrix(data = xt, ncol = 1) %*% t(matrix(data = rt, ncol = 1))
+ # Mise à jour du compteur
+ cc <- cc + 1
+
+ ################################################################################
+ # Dictionary Update Problem : Résolution par block-coordinate descent
+ # with warm restarts
+ # Condition pour la mise à jour selon le choix de la taille du batch
+ if(cc == Tb) {
+ # Boucle pour le choix du nombre d'itération pour la convergence
+ for(i in 1:Tc) {
+ # Mise à jour de la j-ème colonne de Dt "with warm restarts"
+ for(j in 1:k) {
+ # Si At[j,j] == 0, alors Bt[,j] = 0 et At[,j] = 0 par construction
+ # Donc la mise à jour est inutile
+ if(At[j,j]!=0) {
+ # Mise à jour par Gradient Descent
+ uj <- (Bt[,j]-Dt %*% At[,j])/At[j,j] + Dt[,j]
+ # Normalisation pour que Dt appartienne à C (respect de la contrainte sur Dt)
+ Dt[,j] <- uj / max(1,sqrt(sum(uj^2)))
+ }
+ }
+ }
+
+ # Contrôle sur les coefficients de Dt pour éviter les problèmes rencontrés avec le LARS
+ for(i in 1:p) {
+ for(j in 1:k) {
+ Dt[i,j] <- ifelse(Dt[i,j]<10e5*.Machine$double.eps,0,Dt[i,j])
+ }
+ }
+ # Réinitialisation du compteur
+ cc <- 0
+ }
+
+ # Visualisation de l'avancement de l'apprentissage
+ if(Tt>=100 & t %% floor(Tt/100) == 0) cat (".") # cat("Etape t =",t,"\n")
+ if(Tt>=10 & t %% floor(Tt/10) == 0) cat ("|")
+ }
+ cat("\n")
+
+ ################################################################################
+ # Return : Le dictionnaire Dt permettant de bien représenter les features de X
+ res <- list()
+ res$D <- Dt
+ res$D_visual <- plot_mnist_choose(t(Dt),floor(sqrt(k)),
+ title="Atoms Dictionary",
+ subtitle="Visualization of atoms of final dictionary")
+ res$fn <- cost_function(X[ind_ti,], Dt, lambda)
+ res
+}
+
+
+# ----------------------------------------------------------
+# ----------------------------------------------------------
+# visual_structure := Visual Classification
+#
+# Require :
+# X := test dataset (matrix : size = n*p)
+# l := regularization parameter (double : > 0)
+# D := final dictionary (matrix : size = p*k)
+# Y := label test dataset (array : size = n)
+# Tt := update iteration number (int : > 0)
+# nb_rp := nombre de représentants pour une classe (int : > 0)
+# ----------------------------------------------------------
+# ----------------------------------------------------------
+
+visual_structure <- function(X,Y,D,l,Tt,nb_rp) {
+
+ ################################################################################
+ # Initialisation des paramètres n, p, k
+
+ # n := nombre de features de X
+ # p := nombre de parameters de X
+ dim_X <- dim(X) ; n <- dim_X[1] ; p <- dim_X[2]
+ # k := nombre d'atoms de D
+ # p := nombre de parameters des atoms de D
+ dim_D <- dim(D) ; p2 <- dim_D[1] ; k <- dim(D)[2]
+
+ assert_that(p==p2, msg = "Dimension of X and D not compatible")
+
+ ###############################################################################
+ # Initialisation de la matrice de visualisation (ncol = 10 car 10 chiffres)
+ visual <- matrix(data = 0, nrow = k, ncol = 10)
+ col_count <- array(data = 0, dim = c(10))
+
+ cat("0 0.1 0.2 0.3 0.4 ")
+ cat("0.5 0.6 0.7 0.8 0.9 1\n")
+
+ cat ("|")
+ for(t in 1:Tt) {
+ ################################################################################
+ # Extraction d'une feature xt de X
+ xt <- X[t,]
+ yt <- Y[t]
+
+ ################################################################################
+ # Calcul de la décomposition rt de xt
+ try(model <- lars(D, xt, type = "lasso"), silent = T)
+ if(is.na(model)) {
+ next
+ }
+ ind <- which.min(model$lambda>=l)
+ if(is.na(ind)) ind <- model$lambda[1]
+ rt <- coef(model)[ind,]
+
+ ################################################################################
+ # Intégration de la décomposition rt de xt dans la classe k
+ rt <- ifelse(rt==0,0,1)
+ for(i in 1:k) {
+ visual[i,yt+1] <- visual[i,yt+1] + rt[i]
+ }
+ col_count[yt+1] <- col_count[yt+1] + 1
+
+ # Visualisation de l'avancement de l'apprentissage
+ if(Tt>=100 & t %% floor(Tt/100) == 0) cat (".") # cat("Etape t =",t,"\n")
+ if(Tt>=10 & t %% floor(Tt/10) == 0) cat ("|")
+ }
+ cat("\n")
+
+ for(j in 1:10) {
+ visual[,j] <- visual[,j]/col_count[j]
+ }
+
+ repr <- c()
+ for(j in 1:10) {
+ repr <- cbind(repr,Dt[,which(visual[,j]>=sort(visual[,j],decreasing = T)[nb_rp])[1:nb_rp]])
+ }
+
+ res <- list()
+ res$matrix <- visual
+ res$S_visual <- plot_mnist_choose(t(visual),5,
+ title="Frequency of Atoms Dictionary",
+ subtitle="Visualization of frequency of atoms dictionary used")
+ res$R <- repr
+ res$R_visual <- plot_mnist_choose(t(repr),nb_rp,
+ title="k-th decomposition for number",
+ subtitle="Visualization of k-th decomposition for number 0-9")
+ res
+}
+
+# ----------------------------------------------------------
+# Initialisation de X (train)
+{
+ X <- mnist$train$x
+ dim_X <- dim(X)
+ nx <- dim_X[1] ; p <- dim_X[2]
+ Xy <- mnist$train$y
+}
+
+# ----------------------------------------------------------
+# Initialisation de X (test)
+{
+ Y <- mnist$test$x
+ dim_Y <- dim(Y)
+ ny <- dim_Y[1] ; p <- dim_Y[2]
+ Yy <- mnist$test$y
+}
+
+# ----------------------------------------------------------
+# Initialisation de D
+{
+ k <- 100
+ # k <- 200
+
+ # D <- make_dico_1(p,k)
+ # D <- make_dico_2(k)
+ # D <- make_dico_3(p,k)
+ D <- make_dico_4(10) ; k <- 10*10
+}
+
+# ----------------------------------------------------------
+# Paramétrisation de l'algorithme
+{
+ lambda <- 2
+ Tt <- 1000
+ Tb <- 10
+ Tc <- 10
+}
+
+# ----------------------------------------------------------
+# Run
+{
+ RES <- ODL(X, lambda, D, Tt, Tb, Tc)
+ Dt <- RES$D
+ cat(RES$fn,"\n")
+ p1 <- RES$D_visual
+
+ V <- visual_structure(X,Xy,Dt,lambda,1000,5)
+ p2 <- V$S_visual
+ p3 <- V$R_visual
+
+ multiplot(p1,p2,p3,cols = 2, layout=matrix(c(1,3,2,3),nrow=2,byrow=T))
+}
+
+ {
+ colnames(V$matrix) <- 0:9
+ cor_figures <- cor(V$matrix)
+ corrplot.mixed(cor_figures, upper = "square", order = "FPC")
+}
+
+{
+ fn_train <- c()
+ fn_test <- c()
+ llambda <- c(0.01,0.05,0.1,0.25,0.5,1,2)
+ for(l in llambda) {
+ RES <- ODL(X, l, D, 200, Tb, Tc, F)
+ Dt <- RES$D
+ cat(RES$fn,"\n")
+ fn <- cost_function(Y[1:Tt,],D,l)
+ cat(fn,"\n")
+ fn_train <- c(fn_train,RES$fn)
+ fn_test <- c(fn_test,fn)
+ }
+ df <- data.frame(llambda,fn_train,fn_test)
+ ggplot(df, aes(llambda)) + # basic graphical object
+ geom_line(aes(y=fn_train), colour="red") + # first layer
+ geom_line(aes(y=fn_test ), colour="green") # second layer
+}
+
+# ----------------------------------------------------------
+# Initialisation de l'ensenmble de test unitaire
+mnist_test_unitaire <- make_test()
+
+# ----------------------------------------------------------
+# Test sur un élément aléatoire de Y
+{
+ x <- Y[floor(runif(1,1,ny)),]
+ SC <- lars(Dt, x, type = "lasso")
+ ind <- which.min(SC$lambda>=lambda) # max(which(SC$lambda>=lambda))
+ if(is.na(ind)) ind <- SC$lambda[1]
+ alpha <- coef(SC)[ind,]
+ y <- Dt%*%alpha
+ pt <- plot_mnist_choose(matrix(data=c(x,y),nrow = 2, ncol = p, byrow = T),2, subtitle="")
+ pt
+}
+k-sum(alpha==0) ; (k-sum(alpha==0))/k
diff --git a/code/skeleton.R b/code/skeleton.R
new file mode 100755
index 0000000000000000000000000000000000000000..9142ddc5a2f78fb697c80f4c7882b8d7b94a875b
--- /dev/null
+++ b/code/skeleton.R
@@ -0,0 +1,232 @@
+library(ggplot2)
+library(reshape2)
+#setwd("/Volumes/TCDATA/Projets/MNIST")
+setwd("~/Documents/Cours/PRR/MNIST")
+
+# initialize data directory
+data_dir <- "mnist-data"
+dir.create(data_dir, recursive = TRUE, showWarnings = FALSE)
+
+# download the MNIST data sets, and read them into R
+sources <- list(
+ train = list(
+ x = "https://storage.googleapis.com/cvdf-datasets/mnist/train-images-idx3-ubyte.gz",
+ y = "https://storage.googleapis.com/cvdf-datasets/mnist/train-labels-idx1-ubyte.gz"
+ ),
+
+ test = list(
+ x = "https://storage.googleapis.com/cvdf-datasets/mnist/t10k-images-idx3-ubyte.gz",
+ y = "https://storage.googleapis.com/cvdf-datasets/mnist/t10k-labels-idx1-ubyte.gz"
+ )
+)
+
+# read an MNIST file (encoded in IDX format)
+read_idx <- function(file) {
+
+ # create binary connection to file
+ conn <- gzfile(file, open = "rb")
+ on.exit(close(conn), add = TRUE)
+
+ # read the magic number as sequence of 4 bytes
+ magic <- readBin(conn, what = "raw", n = 4, endian = "big")
+ ndims <- as.integer(magic[[4]])
+
+ # read the dimensions (32-bit integers)
+ dims <- readBin(conn, what = "integer", n = ndims, endian = "big")
+
+ # read the rest in as a raw vector
+ data <- readBin(conn, what = "raw", n = prod(dims), endian = "big")
+
+ # convert to an integer vecto
+ converted <- as.integer(data)
+
+ # return plain vector for 1-dim array
+ if (length(dims) == 1)
+ return(converted)
+
+ # wrap 3D data into matrix
+ matrix(converted, nrow = dims[1], ncol = prod(dims[-1]), byrow = TRUE)
+}
+
+mnist <- rapply(sources, classes = "character", how = "list", function(url) {
+
+ # download + extract the file at the URL
+ target <- file.path(data_dir, basename(url))
+ if (!file.exists(target))
+ download.file(url, target)
+
+ # read the IDX file
+ read_idx(target)
+
+})
+
+# convert training data intensities to 0-1 range
+mnist$train$x <- mnist$train$x / 255
+mnist$test$x <- mnist$test$x / 255
+
+# try plotting the pixel intensities for a random sample of 32 images
+plot_mnist_random <- function(X=mnist$test$x, n=36,
+ title="MNIST Image Data",
+ subtitle="Visualization of a sample of images contained in MNIST data set."
+ ) {
+ m <- dim(X)[2]
+ indices <- sample(nrow(X), size = n)
+ data <- array(X[indices, ], dim = c(n, sqrt(m), sqrt(m)))
+ melted <- melt(data, varnames = c("image", "x", "y"), value.name = "intensity")
+ ggplot(melted, aes(x = x, y = y, fill = intensity)) +
+ geom_tile() +
+ scale_fill_continuous(name = "Pixel Intensity") +
+ scale_y_reverse() +
+ facet_wrap(~ image, nrow = sqrt(n), ncol = sqrt(n)) +
+ theme(
+ strip.background = element_blank(),
+ strip.text.x = element_blank(),
+ panel.spacing = unit(0, "lines"),
+ axis.text = element_blank(),
+ axis.ticks = element_blank()
+ ) +
+ labs(
+ title = title,
+ subtitle = subtitle,
+ x = NULL,
+ y = NULL
+ )
+}
+
+plot_mnist_choose <- function(centers, ncol=5,
+ title="MNIST Image Data",
+ subtitle="Visualization of a sample of images contained in MNIST data set."
+ ) {
+ n <- dim(centers)[1]
+ m <- dim(centers)[2]
+ data <- array(centers, dim = c(n, sqrt(m), sqrt(m)))
+ melted <- melt(data, varnames = c("image", "x", "y"), value.name = "intensity")
+ cat("Affichage en cours ...\n")
+ ggplot(melted, aes(x = x, y = y, fill = intensity)) +
+ geom_tile() +
+ scale_fill_continuous(name = "Pixel Intensity") +
+ scale_y_reverse() +
+ facet_wrap(~ image, nrow = n/ncol, ncol = ncol) +
+ theme(
+ strip.background = element_blank(),
+ strip.text.x = element_blank(),
+ panel.spacing = unit(0, "lines"),
+ axis.text = element_blank(),
+ axis.ticks = element_blank()
+ ) +
+ labs(
+ title = title,
+ subtitle = subtitle,
+ x = NULL,
+ y = NULL
+ )
+}
+
+# plot predictions versus actual for small subset
+plot_mnist_pred <- function(data=mnist$test$x, predictions=mnist$test$y, n=30, format=TRUE,
+ title = "MNIST Image Data",
+ subtitle = "Visualization of a sample of images contained in MNIST data set."
+ ) {
+ indices <- sample(nrow(data), n)
+ if(!format) {
+ classes <- vapply(indices, function(i) {
+ predictions$classes[[i]]
+ }, character(1))
+ }
+ else {
+ classes <- as.character(predictions)[indices]
+ }
+
+ data <- array(data[indices, ], dim = c(n, 28, 28))
+ melted <- melt(data, varnames = c("image", "x", "y"), value.name = "intensity")
+ melted$class <- classes
+
+ image_labels <- setNames(
+ sprintf("Predicted: %s\nActual: %s", classes, mnist$test$y[indices]),
+ 1:n
+ )
+
+ ggplot(melted, aes(x = x, y = y, fill = intensity)) +
+ geom_tile() +
+ scale_y_reverse() +
+ facet_wrap(~ image, ncol = 5, labeller = labeller(image = image_labels)) +
+ theme(
+ panel.spacing = unit(0, "lines"),
+ axis.text = element_blank(),
+ axis.ticks = element_blank()
+ ) +
+ labs(
+ title = title,
+ subtitle = subtitle,
+ x = NULL,
+ y = NULL
+ )
+}
+
+mnist_reduction <- function(x) {
+ n <- sqrt(784)
+ xf <- array(x, dim = c(n, n))
+ mnist_red <- matrix(data = 0, nrow = n/2, ncol = n/2)
+ for(i in 1:(n/2)) {
+ for(j in 1:(n/2)) {
+ mnist_red[i,j] <- (xf[1+2*(i-1),1+2*(j-1)] + xf[1+2*(i-1),1+2*(j-1)+1] +
+ xf[1+2*(i-1)+1,1+2*(j-1)] + xf[1+2*(i-1)+1,1+2*(j-1)+1])/4
+ }
+ }
+ array(mnist_red, dim = c((n/2)^2))
+}
+
+mnist_R <- list()
+mnist_R$train$x <- array(data = 0, dim = c(dim(mnist$train$x)[1],14^2))
+mnist_R$test$x <- array(data = 0, dim = c(dim(mnist$test$x)[1],14^2))
+mnist_R$train$y <- mnist$train$y
+mnist_R$test$y <- mnist$test$y
+
+for(i in 1:dim(mnist$train$x)[1]) {
+ mnist_R$train$x[i,] <- mnist_reduction(mnist$train$x[i,])
+}
+for(i in 1:dim(mnist$test$x)[1]) {
+ mnist_R$test$x[i,] <- mnist_reduction(mnist$test$x[i,])
+}
+
+# Multiple plot function
+#
+# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
+# - cols: Number of columns in layout
+# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
+#
+# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
+# then plot 1 will go in the upper left, 2 will go in the upper right, and
+# 3 will go all the way across the bottom.
+#
+multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
+ library(grid)
+
+ # Make a list from the ... arguments and plotlist
+ plots <- c(list(...), plotlist)
+ numPlots = length(plots)
+
+ # If layout is NULL, then use 'cols' to determine layout
+ if (is.null(layout)) {
+ # Make the panel
+ # ncol: Number of columns of plots
+ # nrow: Number of rows needed, calculated from # of cols
+ layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
+ ncol = cols, nrow = ceiling(numPlots/cols))
+ }
+ if (numPlots==1) {
+ print(plots[[1]])
+ } else {
+ # Set up the page
+ grid.newpage()
+ pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
+
+ # Make each plot, in the correct location
+ for (i in 1:numPlots) {
+ # Get the i,j matrix positions of the regions that contain this subplot
+ matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
+ print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
+ layout.pos.col = matchidx$col))
+ }
+ }
+}
diff --git a/papers/DL.pdf b/papers/DL.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..d696c67ed857f6e3a48c79619f60ae01fbe69345
Binary files /dev/null and b/papers/DL.pdf differ
diff --git a/papers/ODLSC.pdf b/papers/ODLSC.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..a05623cc5fcf035aacde45a83583c4360e66c29c
Binary files /dev/null and b/papers/ODLSC.pdf differ
diff --git a/papers/OLMFSC.pdf b/papers/OLMFSC.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..2b19b3721486d946dda153c14bee67bcfa4ad3e9
Binary files /dev/null and b/papers/OLMFSC.pdf differ
diff --git a/papers/SPCA.pdf b/papers/SPCA.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..af29b232f39e7214d356efb435f5714d389f401a
Binary files /dev/null and b/papers/SPCA.pdf differ
diff --git a/report.pdf b/report.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..a8d715495983f9cd682b839dc5409f94eab28999
Binary files /dev/null and b/report.pdf differ