## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2, dpi = 110)
has_glmnet <- requireNamespace("glmnet", quietly = TRUE)
knitr::opts_chunk$set(eval = has_glmnet)
library(glmfsa)
set.seed(20)

## ----nog, eval = !has_glmnet, echo = FALSE, results = "asis"------------------
# cat("**glmnet is not installed, so the comparisons below were not run.**\n")

## ----setup2-------------------------------------------------------------------
library(glmnet)

make_data <- function(n, p, rho = 0.9, k = 10, task = "regression",
                      noise = 0) {
  z <- matrix(rnorm(n * p), n, p); x <- z
  for (j in 2:p) x[, j] <- rho * x[, j - 1] + sqrt(1 - rho^2) * z[, j]
  ti <- seq_len(k) * 10L
  b <- numeric(p); b[ti] <- 1
  if (task == "regression") {
    y <- as.numeric(x %*% b + rnorm(n))
  } else {
    y <- as.integer(x %*% b > 0)
    if (noise > 0) { fl <- sample(n, round(noise * n))
                     y[fl] <- sample(0:1, length(fl), TRUE) }
  }
  list(x = x, y = y, true = ti)
}

## ----accuracy-----------------------------------------------------------------
compare <- function(d, k = 10) {
  f <- fsa(d$x, d$y, k = k)
  g <- glmnet(d$x, d$y)
  i <- which.min(abs(g$df - k))
  gs <- which(as.matrix(coef(g, s = g$lambda[i]))[-1, 1] != 0)
  c(fsa = length(intersect(f$selected[[1]], d$true)) / k,
    glmnet = length(intersect(gs, d$true)) / k)
}

reps <- t(replicate(10, compare(make_data(300, 2000))))
round(colMeans(reps), 3)

## ----bias---------------------------------------------------------------------
d <- make_data(400, 1000)
f <- fsa(d$x, d$y, k = 10)
g <- glmnet(d$x, d$y)
i <- which.min(abs(g$df - 10))

data.frame(
  truth  = 1,
  fsa    = round(mean(coef(f)[d$true + 1, 1]), 3),
  glmnet = round(mean(as.matrix(coef(g, s = g$lambda[i]))[d$true + 1, 1]), 3))

## ----noise--------------------------------------------------------------------
noisy <- function(rep) {
  set.seed(rep)
  d <- make_data(400, 1000, task = "classification", noise = 0.15)
  g <- glmnet(d$x, d$y, family = "binomial")
  i <- which.min(abs(g$df - 10))
  gs <- which(as.matrix(coef(g, s = g$lambda[i]))[-1, 1] != 0)
  yb <- ifelse(d$y == 1, 1, -1)
  c(`fsa-logistic` = length(intersect(
      fsa(d$x, yb, family = "binomial", loss = "logistic",
          k = 10, s = 1e-4)$selected[[1]], d$true)),
    `fsa-lorenz` = length(intersect(
      fsa(d$x, yb, family = "binomial", loss = "lorenz",
          k = 10, s = 1e-4)$selected[[1]], d$true)),
    `glmnet-lasso` = length(intersect(gs, d$true))) / 10
}
round(colMeans(t(sapply(1:8, noisy))), 3)

## ----timing-------------------------------------------------------------------
tm <- function(p) {
  d <- make_data(1000, p)
  kg <- unique(round(exp(seq(log(100), log(2), length.out = 20))))
  c(p = p,
    fsa = system.time(fsa(d$x, d$y, k = kg))[["elapsed"]],
    glmnet = system.time(glmnet(d$x, d$y, nlambda = 20))[["elapsed"]])
}
round(as.data.frame(t(sapply(c(500, 2000, 8000), tm))), 3)

## ----fast---------------------------------------------------------------------
d <- make_data(1000, 8000)
rbind(
  default = system.time(f1 <- fsa(d$x, d$y, k = 10))[1:3],
  fast    = system.time(f2 <- fsa(d$x, d$y, k = 10, mu = 1000,
                                  niter = 150))[1:3])
c(default = length(intersect(f1$selected[[1]], d$true)),
  fast    = length(intersect(f2$selected[[1]], d$true)))

## ----real, eval = has_glmnet--------------------------------------------------
data(dexter)
dim(dexter$x)
sum(Matrix::colSums(dexter$x != 0) == 0)   # columns with no support at all

## ----realcv, eval = has_glmnet------------------------------------------------
set.seed(11)
fid <- sample(rep(1:5, length.out = nrow(dexter$x)))
cv <- cv.fsa(dexter$x, dexter$y, family = "binomial", loss = "lorenz",
             k = c(50, 20, 10, 5), s = c(0, 1e-3, 1e-2, 1e-1, 1),
             foldid = fid, type.measure = "class")
c(k = cv$k.min, s = cv$s.min)

err <- function(p) mean((p > 0) != (as.character(dexter$yvalid) == "pos"))
g <- glmnet(dexter$x, dexter$y, family = "binomial")
pg <- predict(g, dexter$xvalid, type = "link")
sapply(c(5, 10, 20, 50), function(k) {
  f <- fsa(dexter$x, dexter$y, family = "binomial", loss = "lorenz",
           k = k, s = cv$s.min)
  i <- which.min(abs(g$df - k))
  c(glmfsa = err(drop(predict(f, dexter$xvalid, type = "link"))),
    lasso  = err(pg[, i]))
})