params <-
list(family = "lapis", preset = "homage")

## ----setup-opts, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 8,
  fig.height = 4.6,
  fig.align = "center",
  out.width = "96%",
  dpi = 100
)

benchmark_ready <- all(vapply(
  c("bench", "RSpectra", "irlba"),
  requireNamespace,
  logical(1),
  quietly = TRUE
))
benchmark_iterations <- 3L
options(knitr.kable.NA = "not run")

## ----albers-classes, echo=FALSE, results='asis'-------------------------------
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))

## ----helpers, include = FALSE-------------------------------------------------
`%||%` <- function(x, y) if (is.null(x)) y else x

fro_norm <- function(A) {
  if (inherits(A, "sparseMatrix")) {
    sqrt(sum(A@x^2))
  } else {
    sqrt(sum(abs(A)^2))
  }
}

relative_error <- function(x, truth) {
  x <- sort(Re(x), decreasing = TRUE)
  truth <- sort(Re(truth), decreasing = TRUE)
  max(abs(x - truth) / pmax(1, abs(truth)))
}

eigen_backward_error <- function(A, values, vectors) {
  A_dense <- as.matrix(A)
  values <- Re(values)
  vectors <- as.matrix(vectors)
  residual <- A_dense %*% vectors - vectors %*% diag(values, nrow = length(values))
  scale <- fro_norm(A_dense) + abs(values)
  max(sqrt(colSums(abs(residual)^2)) / pmax(scale, .Machine$double.eps))
}

svd_backward_error <- function(A, d, u, v) {
  A_dense <- as.matrix(A)
  d <- Re(d)
  u <- as.matrix(u)
  v <- as.matrix(v)
  left_residual <- A_dense %*% v - u %*% diag(d, nrow = length(d))
  right_residual <- t(A_dense) %*% u - v %*% diag(d, nrow = length(d))
  residual <- sqrt(colSums(abs(left_residual)^2) + colSums(abs(right_residual)^2))
  max(residual / pmax(fro_norm(A_dense) + d, .Machine$double.eps))
}

bench_eval <- function(expr, iterations = benchmark_iterations) {
  expr <- substitute(expr)
  env <- parent.frame()
  last_result <- NULL
  mark <- tryCatch(
    bench::mark(
      last_result <- eval(expr, env),
      iterations = iterations,
      check = FALSE,
      memory = TRUE,
      filter_gc = FALSE
    ),
    error = function(e) e
  )
  if (inherits(mark, "error")) {
    return(list(
      result = NULL,
      median_ms = NA_real_,
      mem_mb = NA_real_,
      error = conditionMessage(mark)
    ))
  }
  list(
    result = last_result,
    median_ms = as.numeric(mark$median[[1L]]) * 1000,
    mem_mb = as.numeric(mark$mem_alloc[[1L]]) / 1024^2,
    error = NA_character_
  )
}

make_dense_hermitian <- function(n, seed = 1L) {
  set.seed(seed)
  X <- matrix(rnorm(n * n), n, n)
  crossprod(X) / n + diag(seq(1, 1.2, length.out = n))
}

make_dense_low_rank <- function(m, n, rank = 8L, noise = 1e-3, seed = 1L) {
  set.seed(seed)
  U <- qr.Q(qr(matrix(rnorm(m * rank), m, rank)))
  V <- qr.Q(qr(matrix(rnorm(n * rank), n, rank)))
  signal <- U %*% diag(seq(rank, 1, length.out = rank), nrow = rank) %*% t(V)
  signal + noise * matrix(rnorm(m * n), m, n)
}

path_laplacian <- function(n) {
  Matrix::bandSparse(
    n,
    k = c(-1L, 0L, 1L),
    diagonals = list(rep(-1, n - 1L), c(1, rep(2, n - 2L), 1), rep(-1, n - 1L))
  )
}

eigen_rows <- function(name, A, k = 6L, tol = 1e-8, seed = 1L,
                       iterations = benchmark_iterations) {
  truth <- eigen(as.matrix(A), symmetric = TRUE, only.values = TRUE)$values[seq_len(k)]
  methods <- c("eigencore", "RSpectra", "base")
  rows <- lapply(methods, function(method) {
    timed <- switch(
      method,
      eigencore = bench_eval({
        set.seed(seed)
        eig_partial(A, k = k, target = largest(), tol = tol)
      }, iterations = iterations),
      RSpectra = bench_eval({
        RSpectra::eigs_sym(A, k = k, which = "LA", opts = list(tol = tol, maxitr = 1000L))
      }, iterations = iterations),
      base = bench_eval({
        eigen(as.matrix(A), symmetric = TRUE)
      }, iterations = iterations)
    )
    if (!is.na(timed$error)) {
      return(data.frame(
        regime = name,
        task = "eigen",
        method = method,
        median_ms = timed$median_ms,
        mem_mb = timed$mem_mb,
        rel_error = NA_real_,
        backward_error = NA_real_,
        residual_check = FALSE,
        eigencore_label = NA_character_,
        status = timed$error,
        stringsAsFactors = FALSE
      ))
    }
    result <- timed$result
    extracted <- switch(
      method,
      eigencore = list(values = values(result), vectors = vectors(result)),
      RSpectra = list(values = result$values, vectors = result$vectors),
      base = {
        ord <- order(result$values, decreasing = TRUE)[seq_len(k)]
        list(values = result$values[ord], vectors = result$vectors[, ord, drop = FALSE])
      }
    )
    backward <- eigen_backward_error(A, extracted$values, extracted$vectors)
    data.frame(
      regime = name,
      task = "eigen",
      method = method,
      median_ms = timed$median_ms,
      mem_mb = timed$mem_mb,
      rel_error = relative_error(extracted$values, truth),
      backward_error = backward,
      residual_check = isTRUE(backward <= tol),
      eigencore_label = if (method == "eigencore") result$method else NA_character_,
      status = "ok",
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

svd_rows <- function(name, A, rank = 6L, tol = 1e-8, seed = 1L,
                     iterations = benchmark_iterations) {
  truth <- svd(as.matrix(A), nu = 0, nv = 0)$d[seq_len(rank)]
  methods <- c("eigencore", "RSpectra", "irlba", "base")
  rows <- lapply(methods, function(method) {
    timed <- switch(
      method,
      eigencore = bench_eval({
        set.seed(seed)
        svd_partial(A, rank = rank, target = largest(), tol = tol)
      }, iterations = iterations),
      RSpectra = bench_eval({
        RSpectra::svds(A, k = rank, nu = rank, nv = rank, opts = list(tol = tol, maxitr = 1000L))
      }, iterations = iterations),
      irlba = bench_eval({
        set.seed(seed)
        irlba::irlba(A, nv = rank, nu = rank, tol = tol)
      }, iterations = iterations),
      base = bench_eval({
        svd(as.matrix(A), nu = rank, nv = rank)
      }, iterations = iterations)
    )
    if (!is.na(timed$error)) {
      return(data.frame(
        regime = name,
        task = "SVD",
        method = method,
        median_ms = timed$median_ms,
        mem_mb = timed$mem_mb,
        rel_error = NA_real_,
        backward_error = NA_real_,
        residual_check = FALSE,
        eigencore_label = NA_character_,
        status = timed$error,
        stringsAsFactors = FALSE
      ))
    }
    result <- timed$result
    extracted <- switch(
      method,
      eigencore = list(d = values(result), u = left_vectors(result), v = right_vectors(result)),
      RSpectra = list(d = result$d, u = result$u, v = result$v),
      irlba = list(d = result$d, u = result$u, v = result$v),
      base = list(d = result$d[seq_len(rank)], u = result$u, v = result$v)
    )
    backward <- svd_backward_error(A, extracted$d, extracted$u, extracted$v)
    data.frame(
      regime = name,
      task = "SVD",
      method = method,
      median_ms = timed$median_ms,
      mem_mb = timed$mem_mb,
      rel_error = relative_error(extracted$d, truth),
      backward_error = backward,
      residual_check = isTRUE(backward <= tol),
      eigencore_label = if (method == "eigencore") result$method else NA_character_,
      status = "ok",
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

metric_table <- function(rows, metric, digits = 4L) {
  regimes <- unique(rows$regime)
  methods <- c("eigencore", "RSpectra", "irlba", "base")
  out <- data.frame(regime = regimes, check.names = FALSE)

  for (method in methods) {
    out[[method]] <- vapply(regimes, function(regime) {
      keep <- rows$regime == regime & rows$method == method & rows$status == "ok"
      if (!any(keep)) NA_real_ else rows[[metric]][which(keep)[1L]]
    }, numeric(1))
    out[[method]] <- signif(out[[method]], digits)
  }

  out
}

timing_table <- function(rows) {
  out <- metric_table(rows, "median_ms")
  method_columns <- setdiff(names(out), "regime")
  out$lowest_median <- vapply(out$regime, function(regime) {
    current <- rows[
      rows$regime == regime & rows$status == "ok",
      c("method", "median_ms"),
      drop = FALSE
    ]
    current$method[[which.min(current$median_ms)]]
  }, character(1))
  out[, c("regime", "lowest_median", method_columns)]
}

quality_table <- function(rows) {
  regimes <- unique(rows$regime)
  do.call(rbind, lapply(regimes, function(regime) {
    current <- rows[rows$regime == regime, , drop = FALSE]
    ok <- current[current$status == "ok", , drop = FALSE]
    data.frame(
      regime = regime,
      successful_methods = sprintf("%d/%d", nrow(ok), nrow(current)),
      max_relative_error = formatC(max(ok$rel_error), format = "e", digits = 2),
      max_backward_error = formatC(max(ok$backward_error), format = "e", digits = 2),
      residual_checks = if (all(ok$residual_check)) "all pass" else "one or more fail",
      stringsAsFactors = FALSE
    )
  }))
}

## ----setup--------------------------------------------------------------------
library(eigencore)

## ----regimes------------------------------------------------------------------
benchmark_regimes <- data.frame(
  regime = c(
    "dense Hermitian",
    "sparse path Laplacian",
    "dense low-rank SVD",
    "tall sparse SVD",
    "wide sparse SVD"
  ),
  input = c("120 x 120 dense", "300 x 300 dgCMatrix", "180 x 70 dense",
            "320 x 60 dgCMatrix", "60 x 320 dgCMatrix"),
  compared_methods = c("eigencore, RSpectra, base",
                       "eigencore, RSpectra, base",
                       "eigencore, RSpectra, irlba, base",
                       "eigencore, RSpectra, irlba, base",
                       "eigencore, RSpectra, irlba, base")
)
knitr::kable(benchmark_regimes)

## ----build-cases, eval = benchmark_ready, include = FALSE---------------------
set.seed(1001)
bench_cases <- list(
  dense_hermitian = make_dense_hermitian(120L, seed = 1001L),
  sparse_laplacian = path_laplacian(300L),
  dense_low_rank_svd = make_dense_low_rank(180L, 70L, rank = 8L, seed = 1002L),
  tall_sparse_svd = Matrix::rsparsematrix(320L, 60L, density = 0.035),
  wide_sparse_svd = Matrix::rsparsematrix(60L, 320L, density = 0.035)
)

## ----run-benchmarks, eval = benchmark_ready, include = FALSE------------------
benchmark_rows <- rbind(
  eigen_rows("dense Hermitian", bench_cases$dense_hermitian,
             k = 6L, seed = 2001L),
  eigen_rows("sparse path Laplacian", bench_cases$sparse_laplacian,
             k = 6L, seed = 2002L),
  svd_rows("dense low-rank SVD", bench_cases$dense_low_rank_svd,
           rank = 6L, seed = 2003L),
  svd_rows("tall sparse SVD", bench_cases$tall_sparse_svd,
           rank = 6L, seed = 2004L),
  svd_rows("wide sparse SVD", bench_cases$wide_sparse_svd,
           rank = 6L, seed = 2005L)
)
stopifnot(all(is.finite(benchmark_rows$median_ms[benchmark_rows$status == "ok"])))
stopifnot(all(is.finite(benchmark_rows$rel_error[benchmark_rows$status == "ok"])))
stopifnot(all(is.finite(benchmark_rows$backward_error[benchmark_rows$status == "ok"])))

## ----timing-table, eval = benchmark_ready, echo = FALSE-----------------------
timing_rows <- timing_table(benchmark_rows)
knitr::kable(
  timing_rows,
  col.names = c("case", "lowest median", "eigencore", "RSpectra", "irlba", "base R"),
  align = c("l", "l", "r", "r", "r", "r"),
  caption = paste("Median solver-call time in milliseconds from",
                  benchmark_iterations, "iterations per method.")
)

## ----timing-summary, eval = benchmark_ready, echo = FALSE, results = 'asis'----
eigencore_lowest <- sum(timing_rows$lowest_median == "eigencore")
cat(sprintf(
  paste0(
    "In this render, eigencore recorded the lowest median in **%d of %d** ",
    "cases. The result is mixed, and the sub-millisecond rows are especially ",
    "sensitive to setup overhead and run-to-run noise."
  ),
  eigencore_lowest,
  nrow(timing_rows)
))

## ----memory-table, eval = benchmark_ready, echo = FALSE-----------------------
memory_rows <- metric_table(benchmark_rows, "mem_mb")
knitr::kable(
  memory_rows,
  col.names = c("case", "eigencore", "RSpectra", "irlba", "base R"),
  align = c("l", "r", "r", "r", "r"),
  caption = "Allocated memory in megabytes."
)

## ----quality-table, eval = benchmark_ready, echo = FALSE----------------------
quality_rows <- quality_table(benchmark_rows)
knitr::kable(
  quality_rows,
  col.names = c("case", "methods completed", "max relative error",
                "max backward error", "residual checks"),
  align = c("l", "c", "r", "r", "l"),
  caption = "Numerical checks across all methods in each case."
)

## ----planner-table, eval = benchmark_ready, echo = FALSE----------------------
planner_rows <- benchmark_rows[
  benchmark_rows$method == "eigencore",
  c("regime", "eigencore_label", "status"),
  drop = FALSE
]
rownames(planner_rows) <- NULL
knitr::kable(
  planner_rows,
  col.names = c("case", "planner label", "status"),
  align = c("l", "l", "l")
)