## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_fixest <- requireNamespace("fixest", quietly = TRUE)

read_csv <- function(file) {
  paths <- c(
    system.file("benchmarks", file, package = "fastconley"),
    file.path("..", "inst", "benchmarks", file),
    file.path("inst", "benchmarks", file)
  )
  paths <- paths[nzchar(paths) & file.exists(paths)]
  if (!length(paths)) stop("Cannot find benchmark file: ", file)
  read.csv(paths[1], stringsAsFactors = FALSE)
}

bench <- read_csv("fastconley-benchmark-results.csv")
large <- read_csv("fastconley-large-data-results.csv")

fmt_sec <- function(x) sprintf("%.3f", x)
fmt_x <- function(x) sprintf("%.1fx", x)
fmt_e <- function(x) formatC(x, format = "e", digits = 2)
fmt_int <- function(x) format(x, big.mark = ",", scientific = FALSE, trim = TRUE)
fmt_pairs <- function(x) formatC(x, format = "e", digits = 2)

session_path <- system.file("benchmarks", "fastconley-large-data-session.txt",
                            package = "fastconley")
if (session_path == "") {
  session_path <- file.path("..", "inst", "benchmarks",
                            "fastconley-large-data-session.txt")
}
session_lines <- if (file.exists(session_path)) readLines(session_path, warn = FALSE) else ""
line_one <- function(pattern) {
  out <- grep(pattern, session_lines, value = TRUE)
  if (length(out)) out[1] else NA_character_
}
strip_value <- function(pattern, prefix) sub(prefix, "", line_one(pattern))

cpu_model <- strip_value("Model name:", "Model name:\\s+")
cpu_count <- strip_value("^CPU\\(s\\):", "CPU\\(s\\):\\s+")
core_count <- strip_value("Core\\(s\\) per socket:", "Core\\(s\\) per socket:\\s+")
thread_count <- strip_value("Thread\\(s\\) per core:", "Thread\\(s\\) per core:\\s+")
mem_line <- line_one("^Mem:")
r_line <- line_one("^R version")
os_line <- strip_value("^Running under:", "Running under:\\s+")
run_date <- sub("\\..*$", "", unique(large$run_date)[1])

## ----motivation-sim, eval=has_fixest------------------------------------------
library(fastconley)
set.seed(42)

n <- 1500
d <- data.frame(lat = runif(n, 40, 46), lon = runif(n, 5, 11))

# Smooth spatial fields via an exponential-covariance Gaussian process.
rad <- pi / 180
D <- 6371 * acos(pmin(1,
  sin(d$lat * rad) %o% sin(d$lat * rad) +
  (cos(d$lat * rad) %o% cos(d$lat * rad)) * cos(outer(d$lon, d$lon, "-") * rad)))
L <- chol(exp(-D / 150) + diag(1e-8, n))
d$x <- drop(crossprod(L, rnorm(n))) + 0.5 * rnorm(n)
e   <- drop(crossprod(L, rnorm(n))) + 0.5 * rnorm(n)
d$y <- 1 + 0.5 * d$x + e

## ----motivation-fit, eval=has_fixest------------------------------------------
library(fixest)
est <- feols(y ~ x, data = d, demeaned = TRUE)

vcov_conley <- function(x) {
  vcovSpHAC(x, lat = "lat", lon = "lon", dist_cutoff = 150,
            kernel = "bartlett", dist_fn = "spherical",
            ncores = 2, data = d)
}

se_tab <- data.frame(
  `standard errors` = c("IID", "Heteroskedasticity-robust", "Conley (150 km)"),
  `se(x)` = sapply(list(vcov(est, vcov = "iid"),
                        vcov(est, vcov = "hetero"),
                        vcov_conley(est)),
                   function(V) sqrt(diag(V))["x"]),
  check.names = FALSE
)
knitr::kable(se_tab, digits = 4, row.names = FALSE)

## ----quick-start-fixest, eval=FALSE-------------------------------------------
# library(fixest)
# library(fastconley)
# 
# est <- feols(
#   outcome ~ treatment + x1 + x2 | county + year,
#   data = d,
#   demeaned = TRUE
# )
# 
# vcov_fastconley <- function(x) {
#   vcovSpHAC(
#     x,
#     lat = "lat",
#     lon = "lon",
#     dist_cutoff = 100,     # kilometers
#     kernel = "uniform",
#     dist_fn = "spherical",
#     ncores = 8,
#     data = d
#   )
# }
# 
# summary(est, vcov = vcov_fastconley)
# etable(est, vcov = vcov_fastconley)

## ----quick-start-fixest-estimation-vcov, eval=FALSE---------------------------
# est <- feols(
#   outcome ~ treatment + x1 + x2 | county + year,
#   data = d,
#   demeaned = TRUE,
#   vcov = vcov_fastconley
# )

## ----quick-start-fixest-matrix, eval=FALSE------------------------------------
# V_conley <- vcov_fastconley(est)
# summary(est, vcov = V_conley)
# etable(est, vcov = V_conley)

## ----quick-start-lfe, eval=FALSE----------------------------------------------
# library(lfe)
# library(fastconley)
# 
# est <- felm(
#   outcome ~ treatment + x1 + x2 | county + year,
#   data = d,
#   keepCX = TRUE
# )
# 
# V_conley <- vcovSpHAC(
#   est,
#   lat = "lat",
#   lon = "lon",
#   dist_cutoff = 100,
#   kernel = "uniform",
#   dist_fn = "spherical",
#   ncores = 8,
#   data = d
# )
# 
# sqrt(diag(V_conley))

## ----quick-start-panel, eval=FALSE--------------------------------------------
# est <- feols(
#   y ~ treatment + x1 + x2 | unit + year,
#   data = panel_data,
#   demeaned = TRUE
# )
# 
# vcov_shac <- function(x) {
#   vcovSpHAC(
#     x,
#     unit = "unit",
#     time = "year",
#     lat = "lat",
#     lon = "lon",
#     dist_cutoff = 500,
#     lag_cutoff = 1,
#     kernel = "bartlett",
#     dist_fn = "spherical",
#     balanced_pnl = TRUE,
#     ncores = 8,
#     data = panel_data
#   )
# }
# 
# summary(est, vcov = vcov_shac)

## ----cutoff-sensitivity, eval=has_fixest--------------------------------------
cutoffs <- c(50, 100, 150, 300, 600)
se_x <- sapply(cutoffs, function(cut) {
  V <- vcovSpHAC(est, lat = "lat", lon = "lon", dist_cutoff = cut,
                 kernel = "bartlett", dist_fn = "spherical",
                 ncores = 2, data = d)
  sqrt(diag(V))["x"]
})
knitr::kable(data.frame(`cutoff (km)` = cutoffs, `se(x)` = round(se_x, 4),
                        check.names = FALSE), row.names = FALSE)

## ----template-point-data, eval=FALSE------------------------------------------
# vcov_points <- function(x) {
#   vcovSpHAC(
#     x,
#     lat = "lat",
#     lon = "lon",
#     dist_cutoff = 100,
#     kernel = "uniform",
#     dist_fn = "spherical",
#     method = "pairwise",
#     ncores = 8,
#     data = d
#   )
# }

## ----template-raster-data, eval=FALSE-----------------------------------------
# vcov_raster <- function(x) {
#   vcovSpHAC(
#     x,
#     lat = "cell_lat",
#     lon = "cell_lon",
#     dist_cutoff = 250,
#     kernel = "bartlett",
#     dist_fn = "spherical",
#     method = "grid",
#     ncores = 8,
#     data = raster_data
#   )
# }

## ----template-panel-data, eval=FALSE------------------------------------------
# vcov_panel <- function(x) {
#   vcovSpHAC(
#     x,
#     unit = "unit",
#     time = "year",
#     lat = "lat",
#     lon = "lon",
#     dist_cutoff = 500,
#     lag_cutoff = 1,
#     kernel = "bartlett",
#     dist_fn = "spherical",
#     balanced_pnl = TRUE,
#     ncores = 8,
#     data = panel_data
#   )
# }

## ----machine-table, echo=FALSE------------------------------------------------
machine_tab <- data.frame(
  item = c("run date", "CPU", "logical CPUs", "cores per socket",
           "threads per core", "memory", "operating system", "R",
           "fastconley", "fixest", "lfe"),
  value = c(run_date, cpu_model, cpu_count, core_count, thread_count, mem_line,
            os_line, r_line, unique(large$fastconley_version)[1],
            unique(large$fixest_version)[1], unique(large$lfe_version)[1]),
  check.names = FALSE
)
knitr::kable(machine_tab)

## ----large-xsection-table, echo=FALSE-----------------------------------------
lx <- large[large$section == "large_xsection", ]
lx <- lx[order(lx$threads), ]
base <- lx$seconds[lx$threads == min(lx$threads)][1]
lx_tab <- data.frame(
  observations = fmt_int(lx$n_obs),
  regressors = lx$k,
  cutoff_km = lx$cutoff_km,
  threads = lx$threads,
  seconds = fmt_sec(lx$seconds),
  speedup_vs_2_threads = fmt_x(base / lx$seconds),
  dense_matrix_size = lx$implied_dense_weight,
  check.names = FALSE
)
knitr::kable(lx_tab)

## ----fixest-comparison-table, echo=FALSE--------------------------------------
xs <- bench[bench$section == "cross_section", ]
xs_fc <- xs[xs$method == "fastconley", ]
xs_fx <- xs[xs$method == "fixest", ]
xs_w <- merge(xs_fc, xs_fx,
              by = c("benchmark", "n_obs", "k", "cutoff_km", "threads"),
              suffixes = c("_fastconley", "_fixest"))
xs_w <- xs_w[order(xs_w$n_obs, xs_w$cutoff_km, xs_w$threads), ]
xs_tab <- data.frame(
  observations = fmt_int(xs_w$n_obs),
  cutoff_km = xs_w$cutoff_km,
  threads = xs_w$threads,
  fastconley_seconds = fmt_sec(xs_w$seconds_fastconley),
  fixest_seconds = fmt_sec(xs_w$seconds_fixest),
  speedup = fmt_x(xs_w$seconds_fixest / xs_w$seconds_fastconley),
  max_abs_difference = fmt_e(xs_w$max_abs_diff_fastconley),
  check.names = FALSE
)
knitr::kable(xs_tab)

## ----panel-table, echo=FALSE--------------------------------------------------
pn <- large[large$section == "large_panel", ]
pn_tab <- data.frame(
  units = fmt_int(pn$n_unit),
  periods = pn$n_time,
  observations = fmt_int(pn$n_obs),
  regressors = pn$k,
  cutoff_km = pn$cutoff_km,
  lag_cutoff = 1,
  threads = pn$threads,
  seconds = fmt_sec(pn$seconds),
  dense_one_period_size = pn$implied_dense_weight,
  check.names = FALSE
)
knitr::kable(pn_tab)

## ----raster-table, echo=FALSE-------------------------------------------------
rg <- large[large$section == "large_grid", ]
rg <- rg[order(rg$kernel), ]
rg_tab <- data.frame(
  kernel = rg$kernel,
  cells = fmt_int(rg$n_obs),
  cutoff_km = rg$cutoff_km,
  threads = rg$threads,
  seconds = fmt_sec(rg$seconds),
  dense_matrix_size = rg$implied_dense_weight,
  implied_within_cutoff_pairs = fmt_pairs(rg$implied_pairs),
  check.names = FALSE
)
knitr::kable(rg_tab)

## ----pixel-table, echo=FALSE--------------------------------------------------
pix <- bench[bench$section == "pixel", ]
pix <- pix[order(pix$cases, decreasing = TRUE), ]
pixel0 <- pix$seconds[pix$method == "pixel=0"][1]
pix_tab <- data.frame(
  setting = pix$method,
  raw_rows = fmt_int(pix$raw_rows),
  coordinate_cases = fmt_int(pix$cases),
  seconds = fmt_sec(pix$seconds),
  speed_vs_pixel0 = fmt_x(pixel0 / pix$seconds),
  max_abs_difference_vs_pixel0 = fmt_e(pix$max_abs_diff),
  check.names = FALSE
)
knitr::kable(pix_tab)

## ----dense-table, echo=FALSE--------------------------------------------------
dense <- bench[bench$section == "dense", ]
dense_fast <- dense[dense$method == "fastconley", ]
dense_base <- dense[dense$method == "dense", ]
dense_fast <- dense_fast[order(dense_fast$n_obs), ]
dense_base <- dense_base[order(dense_base$n_obs), ]
dense_tab <- data.frame(
  observations = fmt_int(dense_fast$n_obs),
  regressors = dense_fast$k,
  cutoff_km = dense_fast$cutoff_km,
  dense_seconds = fmt_sec(dense_base$seconds),
  fastconley_seconds = fmt_sec(dense_fast$seconds),
  speedup = fmt_x(dense_base$seconds / dense_fast$seconds),
  max_abs_difference = fmt_e(dense_fast$max_abs_diff),
  check.names = FALSE
)
knitr::kable(dense_tab)