## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE,
  fig.width = 6,
  fig.height = 4,
  out.width = "100%"
)
library(fastfocal)
library(terra)
library(dplyr)

## -----------------------------------------------------------------------------
library(fastfocal)
library(terra)
library(dplyr)

raster_sizes <- c(100, 250, 500, 1000, 2500, 5000)
kernel_sizes <- seq(100, 1000, 100)
replicates <- 1
res_m <- 30
crs_m <- "EPSG:3857"
set.seed(888)

## -----------------------------------------------------------------------------
rasters <- lapply(raster_sizes, function(size) {
  ext_x <- size * res_m
  ext_y <- size * res_m
  r <- rast(nrows = size, ncols = size, extent = ext(0, ext_x, 0, ext_y), crs = crs_m)
  values(r) <- runif(ncell(r))
  r
})
names(rasters) <- as.character(raster_sizes)

## -----------------------------------------------------------------------------
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 3), mar = c(2, 2, 3, 1))
raster_labels <- paste0(raster_sizes, "x", raster_sizes)
for (i in seq_along(rasters)) {
  plot(rasters[[i]], main = paste("Raster:", raster_labels[i]))
}
par(oldpar)

## ----eval=FALSE---------------------------------------------------------------
# grid <- expand.grid(
#   raster_size = raster_sizes,
#   d = kernel_sizes,
#   method = c("fastfocal", "terra"),
#   stringsAsFactors = FALSE
# )
# 
# dir.create("benchmark_chunks", showWarnings = FALSE)
# 
# benchmark_row <- function(idx) {
#   size <- grid$raster_size[idx]
#   d <- grid$d[idx]
#   method <- grid$method[idx]
#   fname <- sprintf("benchmark_chunks/%s_%d_%dm.csv", method, size, d)
#   if (file.exists(fname)) return(NULL)
#   r <- rasters[[as.character(size)]]
#   times <- sapply(seq_len(replicates), function(i) {
#     t0 <- Sys.time()
#     if (method == "fastfocal") {
#       fastfocal(x = r, d = d, w = "circle", fun = "mean", engine = "auto", pad = "auto")
#     } else {
#       w <- focalMat(r, d, type = "circle")
#       if (all(w == 0)) return(NA_real_)
#       focal(r, w = w, fun = mean, na.rm = TRUE, na.policy = "omit")
#     }
#     as.numeric(difftime(Sys.time(), t0, units = "secs"))
#   })
#   chunk_df <- data.frame(method = method, raster_size = size, d = d, time = times)
#   write.csv(chunk_df, file = fname, row.names = FALSE)
# }
# 
# invisible(sapply(seq_len(nrow(grid)), benchmark_row))
# # After running, you can combine chunks into a single CSV under inst/extdata/benchmark.csv

## -----------------------------------------------------------------------------
bench_path <- system.file("extdata", "benchmark.csv", package = "fastfocal")
if (nzchar(bench_path) && file.exists(bench_path)) {
  df <- read.csv(bench_path)
} else {
  # Fallback tiny demo dataset for CRAN if extdata is not installed
  df <- expand.grid(
    method = c("fastfocal", "terra"),
    raster_size = c(250, 500, 1000),
    d = c(100, 300, 500),
    KEEP.OUT.ATTRS = FALSE,
    stringsAsFactors = FALSE
  )
  set.seed(1)
  df$time <- ifelse(df$method == "fastfocal",
                    runif(nrow(df), 0.05, 0.20),
                    runif(nrow(df), 0.08, 0.35))
}
stopifnot(all(c("method","raster_size","d","time") %in% names(df)))

## -----------------------------------------------------------------------------
# --- summary ---
summary_df <- df %>% 
  group_by(method, raster_size, d) %>% 
  summarize( 
    mean_time = mean(time, na.rm = TRUE), 
    se_time = sd(time, na.rm = TRUE) / sqrt(sum(is.finite(time))), 
    .groups = "drop" 
  ) %>% 
  mutate(raster_label = factor( 
    paste0(raster_size, "x", raster_size), 
    levels = paste0(sort(unique(df$raster_size)), "x", sort(unique(df$raster_size))) 
  )) 

oldpar <- par(no.readonly = TRUE)

layout(matrix(1:6, nrow = 2, byrow = TRUE)) 
par(mar = c(4, 4, 3, 1)) 
cols <- c("fastfocal" = "#0072B2", "terra" = "#D55E00") 
raster_labels <- levels(summary_df$raster_label) 

for (label in raster_labels) { 
  subdf <- subset(summary_df, raster_label == label) 
  if (nrow(subdf) == 0) next 
  plot(NA, 
       xlim = range(subdf$d), 
       ylim = range(subdf$mean_time + subdf$se_time, na.rm = TRUE), 
       xlab = "Kernel size (m)", ylab = "Mean time (s)", 
       main = paste("Raster:", label)) 
  methods <- unique(subdf$method) 
  for (m in methods) { 
    data <- subdf[subdf$method == m, ] 
    lines(data$d, data$mean_time, col = cols[m], type = "b", pch = 16) 
    max_time <- max(subdf$mean_time, na.rm = TRUE) 
    min_se <- 0.001 * max_time 
    se <- ifelse(is.na(data$se_time), 0, data$se_time) 
    se_final <- pmax(se, min_se) 
    suppressWarnings(arrows( 
      x0 = data$d,
      y0 = data$mean_time - se_final,
      x1 = data$d, 
      y1 = data$mean_time + se_final, 
      angle = 90, code = 3, length = 0.05, col = cols[m] 
    )) 
  } 
  legend("topleft", legend = methods, col = cols[methods], pch = 16, lty = 1, bty = "n") 
}
layout(1)
par(oldpar)


## -----------------------------------------------------------------------------
test_r <- rasters[["1000"]]
kernel_d <- 500

# fastfocal
r_fast <- fastfocal(test_r, d = kernel_d, w = "circle", fun = "mean", engine = "auto", pad = "auto")

# terra::focal
w <- focalMat(test_r, kernel_d, type = "circle")
r_terra <- focal(test_r, w = w, fun = mean, na.rm = TRUE, na.policy = "omit")

# Differences
r_diff <- abs(r_fast - r_terra)
v_diff <- values(r_diff)

mean_diff <- mean(v_diff, na.rm = TRUE)
max_diff <- max(v_diff, na.rm = TRUE)

cat("Mean difference:", round(mean_diff, 6), "\n")
cat("Max difference :", round(max_diff, 6), "\n")

## -----------------------------------------------------------------------------
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mar = c(2, 2, 3, 2))
plot(test_r,   main = "Original", col = terrain.colors(20))
plot(r_terra,  main = "terra::focal (500 m)", col = terrain.colors(20))
plot(r_fast,   main = "fastfocal (500 m)", col = terrain.colors(20))
plot(r_diff,   main = "Absolute difference", col = hcl.colors(20, "YlOrRd", rev = TRUE))
par(oldpar)

## -----------------------------------------------------------------------------
sessionInfo()

## -----------------------------------------------------------------------------
citation("fastfocal")