## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 3.5
)

## ----setup--------------------------------------------------------------------
library(lazymatrix)

## ----helper-functions---------------------------------------------------------
library(Matrix)
library(bench)
library(dplyr)
library(ggplot2)
#--------------------------------------------------
# Sparse matrix generator
#--------------------------------------------------
create_sparsematrix_col <- function(n, p, sparsity_col) {
  # sparsity_col = fraction of non-zeros PER COLUMN
  # So each column has round(sparsity_col * n) non-zero entries
  n_nonzero_col <- round(sparsity_col * n)

  i_all <- c()
  j_all <- c()

  for (col in 1:p) {
    i_col <- sample(1:n, n_nonzero_col, replace = FALSE)
    i_all <- c(i_all, i_col)
    j_all <- c(j_all, rep(col, n_nonzero_col))
  }

  pairs <- unique(data.frame(i = i_all, j = j_all))
  x     <- rnorm(nrow(pairs))

  Matrix::sparseMatrix(
    i    = pairs$i,
    j    = pairs$j,
    x    = x,
    dims = c(n, p)
  )
}

#--------------------------------------------------
# Benchmark function Computation Time
#--------------------------------------------------
bench_multiply <- function(sparse_matrix, b) {
  A <- scale(sparse_matrix, center = TRUE, scale = TRUE)
  X <- LazyMatrix(sparse_matrix, "sd", "mean")
  bench::mark(
    dense_product = { A %*% b },
    lazy_product  = { X %*% b },
    check         = FALSE,
    min_iterations = 20
  )
}

#--------------------------------------------------
# Benchmark function for Memory Allocation
#--------------------------------------------------

benchmark_memory <- function(M) {
  dense_size <- tryCatch({
    A_dense <- scale(as.matrix(M), center = TRUE, scale = TRUE)
    mem     <- as.numeric(utils::object.size(A_dense)) / 1024^2
    rm(A_dense); gc()
    data.frame(method = "Dense", success = TRUE, mem_MB = mem)
  }, error = function(e) {
    message("Dense error: ", e$message)
    data.frame(method = "Dense", success = FALSE, mem_MB = NA)
  })

  lazy_size <- tryCatch({
    X_lazy <- LazyMatrix(M, "sd", "mean")
    mem    <- as.numeric(utils::object.size(X_lazy)) / 1024^2
    rm(X_lazy); gc()
    data.frame(method = "Lazy", success = TRUE, mem_MB = mem)
  }, error = function(e) {
    message("Lazy error: ", e$message)
    data.frame(method = "Lazy", success = FALSE, mem_MB = NA)
  })

  bind_rows(dense_size, lazy_size)
}

#--------------------------------------------------
# Plot Theme used throughout the thesis for ggplot2
#--------------------------------------------------
thesis_theme <- function(base_size = 12) {
  theme(
    panel.background = element_rect(fill = "white", colour = NA),
    plot.background  = element_rect(fill = "white", colour = NA),
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    panel.border     = element_rect(colour = "black", fill = NA, linewidth = 0.8),
    text             = element_text(size = base_size),
    axis.text        = element_text(size = rel(1.01)),
    axis.text.x      = element_text(size = rel(1.01), angle = 45, hjust = 1)
  )
}

#--------------------------------------------------
# Plot Theme used throughout vignettes for ggplot2
#--------------------------------------------------

vignette_theme <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
  theme(
    # Subtle background strip for facet labels
    strip.background  = element_rect(fill = "#f0f0f0", colour = NA),
    strip.text        = element_text(face = "bold", size = rel(0.95)),

    # Thin, unobtrusive grid
    panel.grid.major  = element_line(colour = "grey88", linewidth = 0.4),
    panel.grid.minor  = element_blank(),

    # Light border just around facet panels
    panel.border      = element_rect(colour = "grey70", fill = NA, linewidth = 0.5),
    panel.spacing     = unit(0.8, "lines"),

    # Legend inside or on top saves horizontal space in HTML
    legend.position   = "top",
    legend.key.width  = unit(1.8, "lines"),

    # Axis: no need to rotate if labels are short
    axis.text.x       = element_text(size = rel(0.9)),
    axis.text.y       = element_text(size = rel(0.9)),
    axis.title        = element_text(size = rel(0.95)),

    plot.margin       = margin(8, 12, 8, 8)
  )
}

## -----------------------------------------------------------------------------
#--------------------------------------------------
# Parameters (smaller for testing)
#--------------------------------------------------
n              <- 10000
sparsity_cols  <- c(0.05, 0.001, 0.0001)
p_values       <- c(50, 100, 200, 500, 1000, 2000)

#--------------------------------------------------
# Run benchmarks
#--------------------------------------------------
results <- list()
counter <- 1

for (sc in sparsity_cols) {
  cat("Running sparsity_col =", sc, "\n")
  for (p in p_values) {
    cat("  p =", p, "\n")

    M  <- create_sparsematrix_col(n = n, p = p, sparsity_col = sc)
    b  <- rnorm(p)
    bm <- bench_multiply(M, b)

    tmp <- bm %>%
      select(expression, min, median, mem_alloc) %>%
      mutate(sparsity_col = sc, n = n, p = p)

    results[[counter]] <- tmp
    counter <- counter + 1
    rm(M); gc()
  }
}

#--------------------------------------------------
# Process results
#--------------------------------------------------
benchmark_results <- bind_rows(results) %>%
  mutate(
    method     = as.character(expression),
    median_sec = as.numeric(median),
    min_sec    = as.numeric(min)
  )


## -----------------------------------------------------------------------------
#--------------------------------------------------
# Plot for computation time
#--------------------------------------------------
ggplot(
  benchmark_results %>%
    mutate(sc_label = factor(
      paste0("alpha = ", sparsity_col),
      levels = paste0("alpha = ", sort(unique(sparsity_col)))
    )),
  aes(x = p, y = median_sec, color = method, group = method)
) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_log10() +
  facet_wrap(~ sc_label, nrow = 1) +
  labs(
    x     = "Number of columns (p)",
    y     = "Median runtime (seconds, log scale)",
    color = "Method"
  ) +
  vignette_theme(13)

## -----------------------------------------------------------------------------
#--------------------------------------------------
# Parameters
#--------------------------------------------------
n              <- 10000
sparsity_cols  <- c(0.05, 0.001, 0.0001)
p_values       <- c(50, 100, 200, 500, 1000, 2000, 5000)

#--------------------------------------------------
# Process results
#--------------------------------------------------
results <- list()
counter <- 1

for (sc in sparsity_cols) {
  cat("Running sparsity_col =", sc, "\n")
  for (p in p_values) {
    cat("  p =", p, "\n")

    M       <- create_sparsematrix_col(n = n, p = p, sparsity_col = sc)
    mem_res <- benchmark_memory(M)
    mem_res <- mem_res %>%
      mutate(sparsity_col = sc, n = n, p = p)

    results[[counter]] <- mem_res
    counter <- counter + 1
    rm(M); gc()
  }
}

memory_results_col <- bind_rows(results) %>%
  mutate(mem_GB = mem_MB / 1024)


## -----------------------------------------------------------------------------
#--------------------------------------------------
# Plot for Memory Allocation
#--------------------------------------------------
ggplot(
  memory_results_col %>%
    mutate(sc_label = factor(
      paste0("alpha = ", sparsity_col),
      levels = paste0("alpha = ", sort(unique(sparsity_col)))
    )),
  aes(x = p, y = mem_GB, color = method, group = method)
) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_log10(
    labels = scales::label_number(suffix = " GB", accuracy = 0.0001)
  ) +
  facet_wrap(~ sc_label, nrow = 1) +
  labs(
    x        = "Number of columns (p)",
    y        = "Allocated memory (GB, log scale)",
    color    = "Method"
  ) +
  vignette_theme(13)