---
title: "Benchmark"
date: "`r Sys.Date()`"
code-annotations: hover
urlcolor: blue
vignette: >
  %\VignetteIndexEntry{Benchmark}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{quarto::html}
execute:
  eval: !expr identical(tolower(Sys.getenv("NOT_CRAN")), "true")
knitr:
  opts_chunk:
    collapse: true
    comment: "#>"
    out.width: "100%"
editor_options: 
  chunk_output_type: console
---

```{r}
#| label: setup
#| include: false

# Limit threads to avoid a CRAN NOTE
Sys.setenv(OMP_THREAD_LIMIT = 2)

options(scipen = 999)

## Packages needed to *build* the vignette (not to run the benchmarks)
library(ggplot2)
library(patchwork)

## ---------------------------------------------------------------------------
## Where the benchmark results live.
## Works both when the working directory is the package root (interactive use)
## and when it is vignettes/ (what happens at build time).
## ---------------------------------------------------------------------------
bench_dir <- if (dir.exists("vignettes")) "vignettes/data" else "data"
bench_csv <- file.path(bench_dir, "bench-results.csv")

if (!file.exists(bench_csv)) {
  stop(
    "Benchmark results not found at '", bench_csv, "'.\n",
    "Run the benchmark chunks of this vignette locally and execute the ",
    "'export-results' chunk at the bottom before building."
  )
}

bench_results <- utils::read.csv(bench_csv, stringsAsFactors = FALSE)

## ---------------------------------------------------------------------------
## Helpers to turn the stored numbers into prose
## ---------------------------------------------------------------------------

# Subset one operation
bench_op <- function(op) {
  out <- bench_results[bench_results$op == op, ]
  if (nrow(out) == 0) stop("No benchmark results stored for op = '", op, "'")
  out
}

# One value (median seconds or mem_alloc bytes) for one op / size / package
bench_get <- function(op, n, pkg, metric) {
  row <- bench_results[
    bench_results$op == op & bench_results$n == n & bench_results$pkg == pkg,
  ]
  if (nrow(row) != 1) {
    stop(
      "Expected exactly 1 row for op = '", op, "', n = ", n,
      ", pkg = '", pkg, "'; found ", nrow(row),
      ". The results CSV is probably stale."
    )
  }
  row[[metric]]
}

# sf / duckspatial ratio: > 1 means duckspatial wins
bench_ratio <- function(op, n, metric = c("median", "mem_alloc")) {
  metric <- match.arg(metric)
  bench_get(op, n, "sf", metric) / bench_get(op, n, "duckspatial", metric)
}

fmt_num <- function(x) {
  if (x >= 10) format(round(x)) else format(round(x, 1))
}

# "about 50 times faster than {sf}" / "about 1.2 times slower than {sf}"
bench_speed <- function(op, n) {
  r <- bench_ratio(op, n, "median")
  if (r >= 1) {
    paste0("about ", fmt_num(r), " times faster than {sf}")
  } else {
    paste0("about ", fmt_num(1 / r), " times slower than {sf}")
  }
}

# "allocated 18 times less memory" / "allocated 2 times more memory"
bench_mem <- function(op, n) {
  r <- bench_ratio(op, n, "mem_alloc")
  if (r >= 1) {
    paste0("allocated ", fmt_num(r), " times less memory")
  } else {
    paste0("allocated ", fmt_num(1 / r), " times more memory")
  }
}

## ---------------------------------------------------------------------------
## Plotting helpers (now fed by plain numeric columns, not bench objects)
## ---------------------------------------------------------------------------
ggplot_benchmark <- function(
  data,
  log = FALSE,
  show.legend = TRUE,
  ...
  ) {

    sizes <- sort(unique(data$n))
    data$n_lab <- factor(
      data$n,
      levels = sizes,
      labels = format(sizes, big.mark = ".", scientific = FALSE, trim = TRUE)
    )

    p <- ggplot(
      data = data,
      aes(
        x     = n_lab,
        y     = median,
        color = pkg,
        group = pkg
    )) +
      geom_line(
        linewidth   = 0.7,
        linetype    = "dashed",
        alpha       = 0.6,
        show.legend = show.legend
      ) +
      geom_point(
        aes(size = mem_alloc),
        alpha       = 0.9,
        show.legend = show.legend
      ) +
      scale_size_binned(
        name     = "Memory allocated",
        labels   = scales::label_bytes(),
        n.breaks = 4,
        range    = c(1, 8)
      ) +
      scale_color_brewer(palette = "Set1") +
      labs(
        color = "Package",
        x     = "Data size (rows)",
        y     = "Computation time",
        ...
      ) +
      theme_minimal(base_size = 13) +
      theme(
        panel.grid.minor = element_blank(),
        plot.title       = element_text(face = "bold")
      )

    time_lab <- scales::label_number(suffix = " s", accuracy = 0.01)

    if (log) {
      p + scale_y_log10(labels = time_lab)
    } else {
      p + scale_y_continuous(labels = time_lab)
    }
}

ggplot_assemble <- function(plot1, plot2, fun_name) {

  plot1 +
    plot2 +
    plot_annotation(
      title    = paste("Benchmark Comparison for", fun_name),
      subtitle = "Execution time vs. data size",
      theme    = theme(
        plot.title    = element_text(face = "bold", hjust = .5),
        plot.subtitle = element_text(hjust = .5)
      )
    )

}
```

This vignette benchmarks **{duckspatial}** against **{sf}** across several spatial operations, comparing computation time and memory usage as dataset size grows. We plan to extend it with additional operation types in future releases.

The benchmark code below is **not evaluated** when the vignette is built: it takes too long. Instead, it is run locally once before each release, and the results are stored in `vignettes/data/bench-results.csv`, which is what every number and figure in this document is generated from.

Results shown here were produced on `r unique(bench_results$date)` with:

* `duckspatial`: v`r unique(bench_results$v_duckspatial)`

* `duckdb`: v`r unique(bench_results$v_duckdb)`

* `sf`: v`r unique(bench_results$v_sf)`

## TL;DR

{duckspatial} is substantially faster and allocates far less memory than {sf} in almost all cases, with the advantage becoming more pronounced on larger datasets. The one exception is pairwise distance calculation, where the gap in computation time narrows considerably.

## Prepare data

We are going to test the speed of a bunch of functions using simulated data. We use the `make_points()` function to generate `n` random points within the globe, and 10,000 random rectangles.

```{r}
#| eval: false
#| message: false
#| warning: false
#| code-fold: true
#| code-summary: "Set-up"

# Load necessary packages
library(duckspatial)
library(bench)
library(dplyr)
library(sf)
library(ggplot2)
options(scipen = 999)

# Function to generate random points
make_points <- function(n_points) {
    points_df <- data.frame(
      id = 1:n_points,
      x = runif(n_points, min = -180, max = 180),
      y = runif(n_points, min = -90, max = 90),
      value = rnorm(n_points, mean = 100, sd = 15),
      category = sample(c("A", "B", "C", "D"), n_points, replace = TRUE)
  ) |>
    sf::st_as_sf(coords = c("x", "y"), crs = 4326)
}

# Turn a bench::mark() result into plain, CSV-friendly columns:
#  - median    : seconds   (numeric)
#  - mem_alloc : bytes     (numeric)
bench_tidy <- function(bm, op, n) {
  data.frame(
    op        = op,
    n         = n,
    pkg       = bm$pkg,
    median    = as.numeric(bm$median),
    min       = as.numeric(bm$min),
    mem_alloc = as.numeric(bm$mem_alloc),
    n_itr     = bm$n_itr,
    stringsAsFactors = FALSE
  )
}

# Generate datasets of different sizes
withr::with_seed(27, {
  points_sf_100k <- make_points(1e5)
  points_sf_1mi  <- make_points(1e6)
  points_sf_3mi  <- make_points(3e6)
})

# Generate polygons
# Create large polygon dataset (e.g., administrative regions, zones, etc.)
n_polygons <- 10000
polygons_list <- vector("list", n_polygons)

for(i in 1:n_polygons) {
  # Random center point with buffer from edges
  center_x <- runif(1, min = -170, max = 170)
  center_y <- runif(1, min = -80, max = 80)
  
  # Create simple rectangular polygons to avoid geometry issues
  width <- runif(1, min = 0.5, max = 3)
  height <- runif(1, min = 0.5, max = 3)
  
  # Create rectangle coordinates (must be closed: first point = last point)
  x_coords <- c(
    center_x - width/2,
    center_x + width/2,
    center_x + width/2,
    center_x - width/2,
    center_x - width/2  # Close the polygon
  )
  
  y_coords <- c(
    center_y - height/2,
    center_y - height/2,
    center_y + height/2,
    center_y + height/2,
    center_y - height/2  # Close the polygon
  )
  
  # Create polygon matrix
  coords <- cbind(x_coords, y_coords)
  
  # Create polygon (wrapped in list as required by st_polygon)
  polygons_list[[i]] <- st_polygon(list(coords))
}

polygons_sf <- st_sf(
  poly_id    = 1:n_polygons,
  region     = sample(c("North", "South", "East", "West"), n_polygons, replace = TRUE),
  population = sample(1000:1000000, n_polygons, replace = TRUE),
  geometry   = st_sfc(polygons_list, crs = 4326)
)
```

# Spatial join

The spatial join aims to bring the attributes of a dataset to another dataset based on a spatial predicate. One example would be to have an `x` dataset with points that represent observations of wolves. In a `y` dataset, we could have polygons with attributes describing the geographical location (e.g. the name of the country, the name of the region). So, by using a `ST_Join(x, y, "intersects")`, we would assign the attributes of `y` to `x` according to where the observation of the wolf falls.

This operation can be intensive. In the current version of {duckspatial} we see an improvement in speed and memory usage in bigger datasets, as shown in @fig-st-join.

-   Using 1 million points: {duckspatial} was `r bench_speed("join", 1e6)`, and `r bench_mem("join", 1e6)`.

-   Using 3 million points: {duckspatial} was `r bench_speed("join", 3e6)`, and `r bench_mem("join", 3e6)`.

```{r}
#| eval: false
#| message: false
#| code-fold: true
#| code-summary: "Benchmark code - ddbs_join"

# Helper to run the benchmark
run_join_benchmark <- function(points_sf) {
  temp <- bench::mark(
    iterations  = 3,
    check       = FALSE,
    duckspatial = ddbs_join(points_sf, polygons_sf, join = "within"),
    sf          = st_join(points_sf, polygons_sf, join = st_within)
  )
  temp$pkg <- c("duckspatial", "sf")
  bench_tidy(temp, op = "join", n = nrow(points_sf))
}

# Run the benchmark
df_bench_join <- lapply(
  X   = list(points_sf_100k, points_sf_1mi, points_sf_3mi),
  FUN = run_join_benchmark
) |>
  dplyr::bind_rows()
```

```{r}
#| label: fig-st-join
#| fig-cap: "Benchmark that compares spatial join for duckspatial and sf"
#| echo: false
#| warning: false
#| fig-width: 12
#| fig-height: 6

ggplot_assemble(
  plot1 = ggplot_benchmark(
    bench_op("join"), log = FALSE, show.legend = TRUE,
    subtitle = "A) Normal scale"
  ),
  plot2 = ggplot_benchmark(
    bench_op("join"), log = TRUE, show.legend = FALSE,
    subtitle = "B) Log scale"
  ),
  fun_name = "ddbs_join()"
)
```

# Spatial filter

The spatial filter aims to filter rows of `x` based on a spatial relationship with `y`. For example, let's imagine that we have an `x` dataset with observations of wolves all around the world, and we want to filter only those that are in a specific country, for instance in Spain. If the dataset has this attribute, we can just do it with a simple `dplyr::filter()`, but, if the wolves dataset doesn't include a column with the country, we can use a spatial filter. For that, we need a second dataset `y` with the boundaries of Spain, and by using `ST_Filter(x, y, "intersects")`, we would filter only those wolf observations that intersect with Spain's polygon.

In the current version of {duckspatial} we see an improvement in speed and memory usage in bigger datasets, as shown in @fig-st-filter.

-   Using 1 million points: {duckspatial} was `r bench_speed("filter", 1e6)`, and `r bench_mem("filter", 1e6)`.

-   Using 3 million points: {duckspatial} was `r bench_speed("filter", 3e6)`, and `r bench_mem("filter", 3e6)`.

```{r}
#| eval: false
#| message: false
#| code-fold: true
#| code-summary: "Benchmark code - ddbs_filter"

# Helper to run the benchmark
run_filter_benchmark <- function(points_sf) {
  temp <- bench::mark(
    iterations  = 3,
    check       = FALSE,
    duckspatial = ddbs_filter(points_sf, polygons_sf),
    sf          = st_filter(points_sf, polygons_sf)
  )
  temp$pkg <- c("duckspatial", "sf")
  bench_tidy(temp, op = "filter", n = nrow(points_sf))
}

# Run the benchmark
df_bench_filter <- lapply(
  X   = list(points_sf_100k, points_sf_1mi, points_sf_3mi),
  FUN = run_filter_benchmark
) |>
  dplyr::bind_rows()
```

```{r}
#| label: fig-st-filter
#| fig-cap: "Benchmark that compares spatial filter for duckspatial and sf"
#| echo: false
#| warning: false
#| fig-width: 12
#| fig-height: 6

ggplot_assemble(
  plot1 = ggplot_benchmark(
    bench_op("filter"), log = FALSE, show.legend = TRUE,
    subtitle = "A) Normal scale"
  ),
  plot2 = ggplot_benchmark(
    bench_op("filter"), log = TRUE, show.legend = FALSE,
    subtitle = "B) Log scale"
  ),
  fun_name = "ddbs_filter()"
)
```

# Spatial distances

The `ST_Distance(x, y)` calculates the distance between each observation in `x` against each observation in `y`. The default {duckspatial} mode will return a lazy table with three columns: (id_x) the id of the row in `x`; (id_y) the id of the row in `y`; (distance) the actual distance between those pairs of observations. In the case of `mode = 'sf'`, the result will be a sparse matrix. Note that {duckspatial} will use by default the best distance for the input CRS and geometry type.

-   Using 10,000 pairs of points: {duckspatial} was `r bench_speed("distance", 10000)`, and `r bench_mem("distance", 10000)`.

```{r}
#| eval: false
#| message: false
#| code-fold: true
#| code-summary: "Benchmark code - ddbs_distance"

# Helper to run the benchmark
run_distance_benchmark <- function(n) {

  points_sf <- withr::with_seed(27, make_points(n))

  temp <- bench::mark(
    iterations  = 1,
    check       = FALSE,
    duckspatial = ddbs_distance(points_sf, points_sf),
    sf          = st_distance(points_sf, points_sf)
  )
  temp$pkg <- c("duckspatial", "sf")
  bench_tidy(temp, op = "distance", n = n)
}

df_bench_distance <- lapply(
  X   = c(1000, 5000, 10000),
  FUN = run_distance_benchmark
) |>
  dplyr::bind_rows()
```

```{r}
#| label: fig-st-distance
#| fig-cap: "Benchmark that compares spatial distance for duckspatial and sf"
#| echo: false
#| warning: false
#| fig-width: 12
#| fig-height: 6

ggplot_assemble(
  plot1 = ggplot_benchmark(
    bench_op("distance"), log = FALSE, show.legend = TRUE,
    subtitle = "A) Normal scale"
  ),
  plot2 = ggplot_benchmark(
    bench_op("distance"), log = TRUE, show.legend = FALSE,
    subtitle = "B) Log scale"
  ),
  fun_name = "ddbs_distance()"
)
```

# Dissolving geometries

Dissolving geometries consists in merging/aggregating geometries that share a common attribute into a single geometry per group.

In the current version of {duckspatial} we see an improvement in speed and memory usage in bigger datasets, as shown in @fig-st-dissolve.

-   Using 1 million points: {duckspatial} was `r bench_speed("union_agg", 1e6)`, and `r bench_mem("union_agg", 1e6)`.

-   Using 3 million points: {duckspatial} was `r bench_speed("union_agg", 3e6)`, and `r bench_mem("union_agg", 3e6)`.

```{r}
#| eval: false
#| message: false
#| code-fold: true
#| code-summary: "Benchmark code - ddbs_union_agg"

# Helper to run the benchmark
run_union_benchmark <- function(points_sf) {
  temp <- bench::mark(
    iterations  = 3,
    check       = FALSE,
    duckspatial = ddbs_union_agg(points_sf, by = "category"),
    sf          = points_sf |> 
      group_by(category) |> 
      summarise(geometry = st_union(geometry))
  )
  temp$pkg <- c("duckspatial", "sf")
  bench_tidy(temp, op = "union_agg", n = nrow(points_sf))
}

# Run the benchmark
df_bench_union <- lapply(
  X   = list(points_sf_100k, points_sf_1mi, points_sf_3mi),
  FUN = run_union_benchmark
) |>
  dplyr::bind_rows()
```

```{r}
#| label: fig-st-dissolve
#| fig-cap: "Benchmark that compares dissolving geometries by group for duckspatial and sf"
#| echo: false
#| warning: false
#| fig-width: 12
#| fig-height: 6

ggplot_assemble(
  plot1 = ggplot_benchmark(
    bench_op("union_agg"), log = FALSE, show.legend = TRUE,
    subtitle = "A) Normal scale"
  ),
  plot2 = ggplot_benchmark(
    bench_op("union_agg"), log = TRUE, show.legend = FALSE,
    subtitle = "B) Log scale"
  ),
  fun_name = "ddbs_union_agg()"
)
```

# Geometry predicates

Geometry predicates are boolean functions that test spatial relationships between two geometries (they return true or false). They're fundamental to spatial queries, letting you filter or join data based on how geometries relate to one another in space.
Common examples include: ST_Contains, ST_Within, ST_Touches, ST_Overlaps, ST_Crosses, ST_Disjoint, and ST_Equals.

In the following example, we compare `ST_Intersects(geometry A, geometry B)`, which returns true if the two geometries share any point in common (their interiors or boundaries overlap in any way).

-   Using 1 million points: {duckspatial} was `r bench_speed("intersects", 1e6)`, and `r bench_mem("intersects", 1e6)`.

-   Using 3 million points: {duckspatial} was `r bench_speed("intersects", 3e6)`, and `r bench_mem("intersects", 3e6)`.

```{r}
#| eval: false
#| message: false
#| code-fold: true
#| code-summary: "Benchmark code - ddbs_intersects"

# Helper to run the benchmark
run_predicate_benchmark <- function(points_sf) {
  temp <- bench::mark(
    iterations  = 1,
    check       = FALSE,
    duckspatial = ddbs_intersects(points_sf, polygons_sf),
    sf          = st_intersects(points_sf, polygons_sf)
  )
  temp$pkg <- c("duckspatial", "sf")
  bench_tidy(temp, op = "intersects", n = nrow(points_sf))
}

# Run the benchmark
df_bench_predicate <- lapply(
  X   = list(points_sf_100k, points_sf_1mi, points_sf_3mi),
  FUN = run_predicate_benchmark
) |>
  dplyr::bind_rows()
```

```{r}
#| label: fig-st-intersects
#| fig-cap: "Benchmark that compares the geometry predicate ST_Intersects(x, y) for duckspatial and sf"
#| echo: false
#| warning: false
#| fig-width: 12
#| fig-height: 6

ggplot_assemble(
  plot1 = ggplot_benchmark(
    bench_op("intersects"), log = FALSE, show.legend = TRUE,
    subtitle = "A) Normal scale"
  ),
  plot2 = ggplot_benchmark(
    bench_op("intersects"), log = TRUE, show.legend = FALSE,
    subtitle = "B) Log scale"
  ),
  fun_name = "ddbs_intersects()"
)
```

# Export results

Run this chunk once, after all the benchmark chunks above, to refresh the CSV
that this vignette is built from. Everything else in the document (numbers,
prose, figures) updates automatically from it.

```{r}
#| label: export-results
#| eval: false
#| code-fold: true
#| code-summary: "Export the results to CSV"

bench_export <- dplyr::bind_rows(
  df_bench_join,
  df_bench_filter,
  df_bench_distance,
  df_bench_union,
  df_bench_predicate
)

bench_export$date          <- as.character(Sys.Date())
bench_export$v_duckspatial <- as.character(utils::packageVersion("duckspatial"))
bench_export$v_duckdb      <- as.character(utils::packageVersion("duckdb"))
bench_export$v_sf          <- as.character(utils::packageVersion("sf"))

dir.create(bench_dir, recursive = TRUE, showWarnings = FALSE)
utils::write.csv(bench_export, bench_csv, row.names = FALSE)
```