---
title: "Benchmarking fastfocal"
subtitle: "Version 0.1.3 (2025)"
author: "Ho Yi Wan"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Benchmarking fastfocal}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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)
```

## Is it faster?

We compare `fastfocal()` against `terra::focal()` across a range of raster sizes and kernel radii. To keep this vignette fast for CRAN, we load precomputed results if available and provide optional code (disabled) to reproduce full benchmarks locally.

---

## Load libraries and parameters

```{r}
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)
```

---

## Create test rasters

Each raster is square with 30 m resolution.

```{r}
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)
```

---

## Quick peek

```{r}
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)
```

---

## Optional full benchmark (disabled for speed)

Note: running the full grid below can take a while. It is disabled to keep CRAN checks fast. Uncomment to run locally.

```{r 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
```

---

## Load precomputed results (with fallback)

```{r}
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)))
```

---

## Summarize and visualize

```{r}
# --- 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)

```

---

## Bonus: accuracy check

Compare a single case at moderate size.

```{r}
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")
```

### Visual comparison

```{r}
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)
```

---

## Session info

```{r}
sessionInfo()
```

## Citation

To cite the package:

Wan, H. Y. (2025). *fastfocal: Fast Multi-scale Raster Extraction and Moving Window Analysis with FFT*. 
R package version 0.1.3. Zenodo. https://doi.org/10.5281/zenodo.17074691

```{r}
citation("fastfocal")
```