### R code from vignette source 'spareg.Rnw'

###################################################
### code chunk number 1: r setup
###################################################
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
library("ggplot2")
library("spareg")
is_paper <- FALSE
tmpdir <- tempdir()


###################################################
### code chunk number 2: condition_vignette
###################################################
if (is_paper) {
  cat("\\Positivetrue")
} else {
  cat("\\Positivefalse")
}


###################################################
### code chunk number 3: spareg.Rnw:687-689 (eval = FALSE)
###################################################
## install.packages("spareg")
## library("spareg")


###################################################
### code chunk number 4: spareg.Rnw:691-692
###################################################
library("spareg")


###################################################
### code chunk number 5: spareg.Rnw:706-709 (eval = FALSE)
###################################################
## set.seed(1234)
## example_data <- simulate_spareg_data(n = 200, p = 2000, ntest = 200)
## str(example_data)


###################################################
### code chunk number 6: spareg.Rnw:711-714
###################################################
set.seed(1234)
example_data <- simulate_spareg_data(n = 200, p = 2000, ntest = 200)
str(example_data, list.len = 4)


###################################################
### code chunk number 7: first_example
###################################################
set.seed(12)
spar_res <- spar(example_data$x, example_data$y,
  xval = example_data$xtest, yval = example_data$ytest,
  nummods = c(5, 10, 15, 20, 25, 30))
spar_res


###################################################
### code chunk number 8: first_example_cv
###################################################
set.seed(12)
spar_cv <- spar.cv(example_data$x, example_data$y, nfolds = 3,
  nummods = c(5, 10, 15, 20, 25, 30))


###################################################
### code chunk number 9: spareg.Rnw:1239-1240 (eval = FALSE)
###################################################
## coef(spar_res)


###################################################
### code chunk number 10: spareg.Rnw:1242-1243
###################################################
print(coef(spar_res), show = 5L)


###################################################
### code chunk number 11: spareg.Rnw:1245-1246
###################################################
get_intercept(coef(spar_res, nummod = 5, aggregate = "none"))


###################################################
### code chunk number 12: plotmethod
###################################################
library(ggplot2)
p1 <- plot(spar_res)+
  theme(axis.text.x = element_text(angle = 45,vjust = 0.5))
p2 <- plot(spar_res, plot_type = "val_numactive") +
  theme(axis.text.x = element_text(angle = 45,vjust = 0.5))
p4 <- plot(spar_res, plot_type = "res_vs_fitted", xfit = example_data$xtest,
  yfit = example_data$ytest)
p3 <- plot(spar_res,  plot_type = "coefs")
p <- ggpubr::ggarrange(p1, p2, p3, p4,
  ncol = 2, nrow = 2)
p


###################################################
### code chunk number 13: parallel_example (eval = FALSE)
###################################################
## ex4 <- simulate_spareg_data(n = 1000, p = 2000, ntest = 1000, seed = 123)
## library(doParallel)
## library(doRNG)
## cl <- makeCluster(2, type = "PSOCK")
## registerDoParallel(cl)
## registerDoRNG(seed = 123)
## spar_res_par <- spar(ex4$x, ex4$y, screencoef = screen_cor(),
##   rp = rp_gaussian(), nummods = 50, parallel = TRUE)
## stopCluster(cl)


###################################################
### code chunk number 14: spareg.Rnw:1407-1414
###################################################
generate_scr_sirs <- function(object, x, y, ...) {
  ctrl <- object$control[names(object$control) %in%
                           names(formals(VariableScreening::screenIID))]
  res_screen <- do.call(function(...)
    VariableScreening::screenIID(x, y, ...), ctrl)
  res_screen$measurement
}


###################################################
### code chunk number 15: sirs_chunk
###################################################
screen_sirs <- constructor_screencoef(generate_fun = generate_scr_sirs,
  name = "screen_sirs")


###################################################
### code chunk number 16: screen_sirs
###################################################
set.seed(123)
spar_example <- spar(example_data$x, example_data$y,
  screencoef = screen_sirs(type = "fixed", control = list(method = "SIRS")),
  rp = rp_sparse(psi = 1/sqrt(ncol(example_data$x))), measure = "mse")


###################################################
### code chunk number 17: spareg.Rnw:1449-1453
###################################################
simulate_haar <- function(m, p) {
  R0 <- matrix(1/sqrt(p) * rnorm(p * m), nrow = p, ncol = m)
  t(qr.Q(qr(R0), complete = FALSE))
}


###################################################
### code chunk number 18: update_cannings
###################################################
update_rp_cannings <- function(object, x, y, family, ...) {
  if (is.null(object$control$alpha)) object$control$alpha <- 0
  if (is.null(object$control$B))  object$control$B <- 50
  if (is.null(object$control$xi)) object$control$xi <- 0.25
  if (is.null(object$control$family)) object$control$family <- family
  object
}


###################################################
### code chunk number 19: generate_cannings
###################################################
generate_cannings <- function(object, x, y, m, included_vector, ...) {
  xs <- x[, included_vector]; n <- nrow(x);  p <- ncol(xs)
  id_test <- sample(n, size = n * object$control$xi)
  xtrain <- xs[-id_test, ];  xtest <- xs[id_test,]
  ytrain <- y[-id_test];  ytest <- y[id_test]
  control_glm <-
    object$control[names(object$control) %in% names(formals(stats::glm.fit))]
  best_val <- Inf
  for (b in seq_len(object$control$B)) {
    RM <- simulate_haar(m, p)
    xrp <- tcrossprod(xtrain, RM)
    mod <- do.call(function(...)
      glm.fit(x = cbind(1, xrp), y = ytrain, ...), control_glm)
    eta_test <- cbind(1, tcrossprod(xtest, RM)) %*% coef(mod)
    pred <- object$control$family$linkinv(eta_test)
    out_dev <- sum(object$control$family$dev.resid(ytest, pred, wt = 1))
    if (out_dev < best_val) {
      best_val <- out_dev; best_RM <- RM
    }
  }
  return(best_RM)
}


###################################################
### code chunk number 20: define_rp_cannings
###################################################
rp_cannings <- constructor_randomprojection(
  name = "RP Cannings and Samworth 2017",
  generate_fun = generate_cannings, update_fun = update_rp_cannings)


###################################################
### code chunk number 21: data_cannings
###################################################
set.seed(1234)
ex2 <- simulate_spareg_data(n = 200, p = 500, ntest = 100)
ys <- (ex2$y > 0) + 0; ysval <- (ex2$ytest > 0) + 0


###################################################
### code chunk number 22: rp_cannings_1
###################################################
set.seed(1234)
spar_cs <- spar(x = ex2$x, y = ys, family = binomial(),
  rp = rp_cannings(msup = 20), model = spar_glm(),
  nus = 0, nummods = 50, xval = ex2$xtest,  yval = ysval)


###################################################
### code chunk number 23: rp_cannings_2
###################################################
set.seed(1234)
spar_cw <- spar(x = ex2$x, y = ys, family = binomial(),
  rp = rp_cw(data = FALSE, msup = 20), model = spar_glm(),
  nus = 0, nummods = 50, xval = ex2$xtest, yval = ysval)


###################################################
### code chunk number 24: compare
###################################################
rbind("cannings" = get_measure(spar_cs), "cw" = get_measure(spar_cw))


###################################################
### code chunk number 25: spareg.Rnw:1568-1578
###################################################
model_glmrob <- function(object, x, y, z, ...) {
  requireNamespace("robustbase")
  fam <- object$control$family
  f <- ifelse(fam$family == "gaussian" & fam$link == "identity",
    robustbase::lmrob, robustbase::glmrob)
  glmrob_res <- do.call(function(...)
    f(y ~ as.matrix(z), ...), object$control)
  cfs <- coef(glmrob_res)
  list(gammas = cfs[-1], intercept = cfs[1])
}


###################################################
### code chunk number 26: spar_glmrob
###################################################
spar_glmrob <- constructor_sparmodel(generate_fun = model_glmrob)


###################################################
### code chunk number 27: generate_example_data3
###################################################
set.seed(123)
example_data3 <- simulate_spareg_data(n = 100, p = 1000,
  ntest = 100, snr = 10, beta_vals = c(-1, 1)/10)
ypois <- round(exp(example_data3$y)) + 1; x <- example_data3$x;
np <- ncol(x) * nrow(x)
id_outliers_x <- sample(seq_len(np),  0.25 * np)
x[id_outliers_x] <- x[id_outliers_x] + 50


###################################################
### code chunk number 28: run_glmrob_glm
###################################################
set.seed(1234)
spar_rob_res <- spar(x, ypois, family = poisson(),
  model = spar_glmrob(), rp = rp_gaussian(msup = 25), measure = "mae")
set.seed(1234)
spar_res <- spar(x, ypois, family = poisson(),
  model = spar_glm(), rp = rp_gaussian(msup = 25), measure = "mae")


###################################################
### code chunk number 29: compare_glmrob_glm
###################################################
best_rob <- get_model(spar_rob_res, opt_par = "best")
best_glm <- get_model(spar_res, opt_par = "best")


###################################################
### code chunk number 30: compare_glmrob_glm
###################################################
rbind("robglm" = get_measure(best_rob), "glm" = get_measure(best_glm))


###################################################
### code chunk number 31: urls
###################################################
if (is_paper) {
  url1 <- "https://web.archive.org/web/20150922051706/"
  url2 <- "http://isomap.stanford.edu/face_data.mat.Z"
  ## On Ubuntu or MacOS use code below
  ## On Windows, download locally and unzip
  if (!file.exists(file.path(tmpdir, "face_data.mat"))) {
    # system('uncompress face_data.mat.Z')
    library("R.matlab")
    download.file(paste0(url1, url2),
                  file.path(tmpdir, "face_data.mat.Z"))
    system(sprintf('uncompress %s', file.path(tmpdir, "face_data.mat.Z")))
  }
}


###################################################
### code chunk number 32: unzip_file
###################################################
if (is_paper) {
  library("R.matlab")
  facedata <- readMat(file.path(tmpdir,"face_data.mat"))
  x <- t(facedata$images)
  y <- facedata$poses[1,]
  x[, apply(x, 2, sd) < 0.01] <- 0
}


###################################################
### code chunk number 33: read_matlab (eval = FALSE)
###################################################
## library("R.matlab")
## facedata <- readMat("face_data.mat")
## x <- t(facedata$images); y <- facedata$poses[1, ]
## x[, apply(x, 2, sd) < 0.01] <- 0


###################################################
### code chunk number 34: std_x
###################################################
if (is_paper) {
  set.seed(123)
  ntot <- length(y); ntest <- ntot * 0.25
  testind <- sample(ntot, ntest, replace = FALSE)
  xtrain <- as.matrix(x[-testind, ]); ytrain <- y[-testind]
  xtest <- as.matrix(x[testind, ]); ytest <- y[testind]
}


###################################################
### code chunk number 35: spareg.Rnw:1693-1697 (eval = FALSE)
###################################################
## set.seed(123)
## testind <- sample(length(y), length(y) * 0.25, replace = FALSE)
## xtrain <- as.matrix(x[-testind, ]); ytrain <- y[-testind]
## xtest <- as.matrix(x[testind, ]); ytest <- y[testind]


###################################################
### code chunk number 36: spareg.Rnw:1703-1724
###################################################
if (is_paper) {
  file <- file.path(tmpdir, "faces_res.rda")
  if (!file.exists(file)) {
    set.seed(123)
    control_glmnet <- list(lambda.min.ratio = 0.001)
    library("spareg")
    Sys.time()
    spar_faces <- spar.cv(
      xtrain, ytrain,
      model = spar_glm(),
      screencoef = screen_glmnet(control = control_glmnet,
                                 reuse_in_rp = TRUE),
      rp = rp_cw(data = TRUE),
      nummods = c(5, 10, 20, 50),
      measure = "mse")
    Sys.time()
    save(spar_faces, file = file)
  } else {
    load(file)
  }
}


###################################################
### code chunk number 37: spareg.Rnw:1726-1732 (eval = FALSE)
###################################################
## set.seed(123)
## control_glmnet <- list(lambda.min.ratio = 0.001)
## spar_faces <- spar.cv(xtrain, ytrain,
##   model = spar_glm(), rp = rp_cw(data = TRUE),
##   screencoef = screen_glmnet(control = control_glmnet, reuse_in_rp = TRUE),
##   nummods = c(5, 10, 20, 50), measure = "mse")


###################################################
### code chunk number 38: spareg.Rnw:1740-1741 (eval = FALSE)
###################################################
## spar_faces


###################################################
### code chunk number 39: spareg.Rnw:1743-1744
###################################################
if (is_paper) spar_faces


###################################################
### code chunk number 40: plotfacesmeasure (eval = FALSE)
###################################################
## plot(spar_faces)
## plot(spar_faces, nummod = 5)


###################################################
### code chunk number 41: plotfacesmeasurereal_best
###################################################
if (is_paper) {
  ## best
  p <- plot(spar_faces, digits = 3L) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                     hjust = 1))
  p
}


###################################################
### code chunk number 42: plotfacesmeasurereal_1se
###################################################
if (is_paper) {
  ## 1se
  p2 <- plot(spar_faces, nummod=5, digits = 3L) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                     hjust = 1))
  p2
}


###################################################
### code chunk number 43: facecoef1se (eval = FALSE)
###################################################
## face_coef_1se <- coef(spar_faces, opt_par = "1se")
## summary(face_coef_1se)


###################################################
### code chunk number 44: facecoef1se
###################################################
if (is_paper) {
  face_coef_1se <- coef(spar_faces, opt_par = "1se")
  summary(face_coef_1se)
}


###################################################
### code chunk number 45: tarp_comparison_true (eval = FALSE)
###################################################
## set.seed(123)
## m_l <- 2 * log(ncol(xtrain)); m_u <- 3/4 * nrow(xtrain)
## tarp_faces <-  spar.cv(xtrain, ytrain, screencoef = screen_cor(),
##   rp = rp_sparse(psi = 1/3, mslow = m_l, msup = m_u), model = spar_glm(),
##   nus = 0, nummods = c(5, 10, 20, 50), measure = "mse")


###################################################
### code chunk number 46: tarp_comparison_true
###################################################
if(is_paper) {
  set.seed(123)
  m_l <- 2 * log(ncol(xtrain)); m_u <- 3/4 * nrow(xtrain)
  tarp_faces <-  spar.cv(xtrain, ytrain,
                         model = spar_glm(),
                         screencoef = screen_cor(),
                          nus = 0,
                         rp = rp_sparse(psi = 1/3, mslow = m_l, msup = m_u),
                         nummods = c(5, 10, 20, 50), measure = "mse")
  ynew_tarp <- predict(tarp_faces, xnew = xtest, opt_par = "1se")
}


###################################################
### code chunk number 47: test_tarp_comparison_false (eval = FALSE)
###################################################
## ynew_spar <- predict(spar_faces, opt_par = "1se", xnew = xtest)
## ynew_tarp <- predict(tarp_faces, opt_par = "1se", xnew = xtest)
## c("SPAR" = mean((ytest - ynew_spar)^2),
##   "TARP" = mean((ytest - ynew_tarp)^2))


###################################################
### code chunk number 48: test_tarp_comparison_true
###################################################
if(is_paper) {
ynew_spar <- predict(spar_faces, opt_par = "1se", xnew = xtest)
ynew_tarp <- predict(tarp_faces, opt_par = "1se", xnew = xtest)
  c("SPAR" = mean((ytest - ynew_spar)^2),
    "TARP" = mean((ytest - ynew_tarp)^2))}


###################################################
### code chunk number 49: plotfacesobs179
###################################################
if (is_paper) {
  i <- 179
  p <- ggplot(data.frame(X = rep(1:64,each=64),
                         Y = rep(64:1,64),
                         Z = facedata$images[,i]),
              aes(X, Y, fill = Z)) +
    geom_tile() +
    theme_void() +
    ggtitle(paste0("y = ", round(facedata$poses[1, i],1))) +
    theme(legend.position = "none",
          plot.title = element_text(hjust = 0.5))
  p
}


###################################################
### code chunk number 50: plotobsvspred179
###################################################
if (is_paper) {
  id <- 3
  p4 <- ggplot(data.frame(X = rep(1:64, each = 64),
                          Y = rep(64:1, 64),
                          effect = xtest[id,] * face_coef_1se$beta),
               aes(X, Y, fill = effect)) +
    geom_tile() +
    theme_void() +
    scale_fill_gradient2() +
    ggtitle(bquote(hat(y) == .(round(ynew_spar[id]))), ) +
    theme(plot.title = element_text(hjust = 0.5))
  p4
}


###################################################
### code chunk number 51: readdarwinreal
###################################################
if (is_paper) {
  # if (file.exists(file.path(tmpdir, "data.zip"))) {
  #  unzip(file.path(tmpdir, "data.zip")exdir = )
  #} else {
  if (!file.exists(file.path(tmpdir, "darwin.zip"))) {
    download.file("https://archive.ics.uci.edu/static/public/732/darwin.zip",
                  file.path(tmpdir, "darwin.zip"))
    unzip(zipfile = file.path(tmpdir, "darwin.zip"),
          files = "data.csv",
          exdir = tmpdir)
  }
  darwin_tmp <- read.csv(file.path(tmpdir, "data.csv"),
                         stringsAsFactors = TRUE)
}


###################################################
### code chunk number 52: readdarwinzip (eval = FALSE)
###################################################
## darwin_tmp <- read.csv("data.csv", stringsAsFactors = TRUE)


###################################################
### code chunk number 53: darwinimpute (eval = FALSE)
###################################################
## darwin_orig <- list(y = as.numeric(darwin_tmp$class) - 1,
##   x = darwin_tmp[, !(colnames(darwin_tmp) %in% c("ID", "class"))])
## tmp <- cellWise::DDC(as.matrix(darwin_orig$x),
##   list(returnBigXimp = TRUE, tolProb = 0.999, silent = TRUE))
## darwin <- list(x = tmp$Ximp, y = darwin_orig$y)


###################################################
### code chunk number 54: spareg.Rnw:1933-1942
###################################################
if (is_paper) {
  darwin_orig <- list(
    x = darwin_tmp[, !(colnames(darwin_tmp) %in% c("ID", "class"))],
    y = as.numeric(darwin_tmp$class) - 1)
  tmp <- suppressMessages(cellWise::DDC(
    as.matrix(darwin_orig$x),
    list(returnBigXimp = TRUE, tolProb = 0.999, silent = TRUE)))
  darwin <- list(x = tmp$Ximp, y = darwin_orig$y)
}


###################################################
### code chunk number 55: spardarwin (eval = FALSE)
###################################################
## set.seed(1234)
## spar_darwin <- spar.cv(darwin$x, darwin$y, family = binomial(logit),
##   screencoef = screen_glmnet(reuse_in_rp = TRUE), nummods = c(10, 20, 30, 50),
##   measure = "1-auc")
## spar_darwin


###################################################
### code chunk number 56: spardarwinreal
###################################################
if (is_paper) {
  file <- file.path(tmpdir, "darwin_res.rda")
  if (!file.exists(file)) {
    set.seed(1234)
    spar_darwin <- spareg::spar.cv(darwin$x, darwin$y,
                                   family = binomial(logit),
                                   screencoef = screen_glmnet(reuse_in_rp = TRUE),
                                   nummods = c(10, 20, 30, 50),
                                   measure = "1-auc")
    save(spar_darwin, file = file)
  } else {
    load(file = file)
  }
  spar_darwin
}


###################################################
### code chunk number 57: plotdarwinact (eval = FALSE)
###################################################
## plot(spar_darwin, plot_type = "val_numactive")


###################################################
### code chunk number 58: plotdarwinactnumactive
###################################################
if (is_paper) {
  p <- plot(spar_darwin, plot_type = "val_numactive", digits=3L)
  p
}


###################################################
### code chunk number 59: plotdarwincoef
###################################################
if (is_paper) {
  ntasks <- 25
  nfeat <- 18
  reorder_ind <- c(outer(
    (seq_len(ntasks) - 1) * nfeat,
    seq_len(nfeat), "+"))
  feat_names <- sapply(colnames(darwin$x)[seq_len(nfeat)],
                       function(name) substr(name, 1, nchar(name) - 1))
  p <- plot(spar_darwin,"coefs",coef_order = reorder_ind) +
    geom_vline(xintercept = 0.5 + seq_len(ntasks - 1) * ntasks,
               alpha = 0.2, linetype = 2) +
    annotate("text",x = (seq_len(nfeat) - 1) * ntasks + 12,
             y = 40,
             label=feat_names, angle = 90,
             size = 3)
  p
}


###################################################
### code chunk number 60: session_info
###################################################
sessionInfo()