---
title: "Datashading: millions of points, honestly"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Datashading: millions of points, honestly}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 5,
  fig.height = 5,
  dpi = 100,
  dev = "png"
)
library(vellum)
```

Past a few hundred thousand points, drawing one marker each stops working.
Markers pile on top of markers, the file balloons, and a solid blob hides
whatever structure the data had. Drawing the markers *faster* does not fix this,
because the picture itself is wrong: overplotting throws away density.

`datashade()` takes the approach [datashader](https://datashader.org) made
popular: do not draw markers at all. Bin the points into a canvas-sized grid in
one pass, count how many land in each cell, and colour each cell by that count.
The cost is decoupled from both the number of points and the amount of
overplotting, and the result is honest about density. `datashade()` returns a
single `raster_grob()` you draw like anything else.

The examples here are strange attractors, the same gallery shipped in
`system.file("examples/attractors.R", package = "vellum")`. They are a good
stress test: a chaotic map iterated for millions of points, with density
spanning many orders of magnitude between the sparse outer wisps and the bright
folded core.

## Generating an orbit

An attractor orbit is sequential: each point is a function of the one before it,
so it cannot be vectorised in R. vellum ships a small Rust kernel for it,
`rs_attractor()`, which returns 10 million points in a fraction of a second.
It is internal (hence the `:::`), exposed here only to feed the example; the
subject of this article is what happens *after* the cloud exists.

```{r}
#| label: orbit
n_points <- 2e6

attractor <- function(kind, p, n = n_points, x0 = 0.1, y0 = 0.1) {
  n <- as.integer(n)
  v <- vellum:::rs_attractor(kind, n, p[1], p[2], p[3], p[4], x0, y0)
  list(x = v[seq_len(n)], y = v[n + seq_len(n)])
}

pts <- attractor("clifford", c(-1.4, 1.6, 1.0, 0.7))
length(pts$x)
```

## Shading in one call

Hand the cloud to `datashade()`. `width` and `height` are the aggregation grid
in cells, which is also the output raster in pixels, and `colors` is the
low-to-high density ramp. A white first stop lets the sparse regions fade into
the page.

```{r}
#| label: clifford
blues <- c("#ffffff", "#c6dbef", "#6baed6", "#2171b5", "#08306b")

vl_scene(5, 5, bg = "white") |>
  draw(datashade(
    pts$x, pts$y,
    width = 500, height = 500,
    colors = blues
  ))
```

Two million points, and the work scaled with the half-million grid cells, not
the points. The same call handles twenty million with no change to the picture's
size on disk: still one 500x500 raster.

## Why aggregation wins

The whole cost is one linear pass over the points to bin them, then a colour
lookup per grid cell. Concretely, `datashade()` calls an aggregation kernel that
is the O(N) heart of the method:

```{r}
#| label: aggregate
counts <- vellum:::rs_aggregate_2d(
  pts$x, pts$y,
  NULL,          # optional per-point weights
  120L, 120L,    # grid size
  min(pts$x), max(pts$x), min(pts$y), max(pts$y)
)
dim(counts)
range(counts)     # the densest cell holds this many points
mean(counts > 0)  # fraction of the grid the orbit actually touches
```

That count grid is all the downstream cost depends on. Whether it came from two
million points or two hundred million, the shading step sees the same 120x120
matrix. Overplotting, the thing that ruins a scatter of markers, is exactly what
the count *measures* here rather than hides.

## Mapping density to colour

The densest cell of an attractor can hold thousands of times more points than a
faint outer cell. Map that range linearly and the bright core saturates while
everything faint collapses to nearly the background. The `how` argument controls
the density-to-colour mapping; datashader's default, `"eq_hist"` (histogram
equalisation), allocates colour by rank so structure stays visible across the
whole range.

```{r}
#| label: how
#| fig.width: 6
#| fig.height: 6
ramp <- c("#ffffff", "#08306b")

panel <- function(how) {
  datashade(pts$x, pts$y, width = 300, height = 300, colors = ramp, how = how)
}

vl_scene(6, 6, bg = "white") |>
  push(vl_viewport(layout = grid_layout(
    widths = vl_unit(c(1, 1), "null"),
    heights = vl_unit(c(1, 1), "null")
  ))) |>
  push(vl_viewport(row = 1, col = 1)) |> draw(panel("linear"))  |> pop() |>
  push(vl_viewport(row = 1, col = 2)) |> draw(panel("log"))     |> pop() |>
  push(vl_viewport(row = 2, col = 1)) |> draw(panel("cbrt"))    |> pop() |>
  push(vl_viewport(row = 2, col = 2)) |> draw(panel("eq_hist")) |> pop()
```

Clockwise from top-left: `"linear"` (only the core survives), `"log"`, `"cbrt"`,
and `"eq_hist"` (the fullest picture). `"log"` and `"cbrt"` are useful middle
grounds when you want a mapping with a fixed analytic meaning rather than one
that depends on the data's rank distribution.

## Lining up with data axes

`datashade()` bins over `xlim` by `ylim` and returns a raster that fills its
viewport (`npc` `0..1`). To place it against data axes, draw it inside a
`vl_viewport()` whose `xscale` / `yscale` match the same limits, then draw axes,
labels, or reference lines in `"native"` units in that viewport. For crisp bins,
match `width` / `height` to the viewport's pixel size and keep the default
`interpolate = FALSE`.

For the attractors we centre a square window on the orbit so a square cell
maps without distortion:

```{r}
#| label: window
window <- function(x, y, pad = 1.05) {
  xr <- range(x)
  yr <- range(y)
  half <- max(diff(xr), diff(yr)) / 2 * pad
  list(
    xlim = mean(xr) + c(-half, half),
    ylim = mean(yr) + c(-half, half)
  )
}

w <- window(pts$x, pts$y)
str(w)
```

Passing that window as `xlim` / `ylim` keeps the aspect honest regardless of how
the raw orbit range differs between the two axes.

## A gallery

Because each shaded attractor is just a grob, a gallery is a layout of them.
Push a `grid_layout()` of `"null"` (flexible) tracks and drop one `datashade()`
raster into each cell. Square cells (from the square window above) keep every
map undistorted.

```{r}
#| label: gallery
#| fig.width: 6
#| fig.height: 6
gallery <- list(
  list(kind = "clifford", p = c(-1.4,  1.6,  1.0,  0.7),
       pal = c("#ffffff", "#c6dbef", "#6baed6", "#2171b5", "#08306b")),
  list(kind = "dejong",   p = c( 1.4, -2.3,  2.4, -2.1),
       pal = c("#ffffff", "#fdd0a2", "#fd8d3c", "#d94801", "#7f2704")),
  list(kind = "svensson", p = c( 1.5, -1.8,  1.6,  0.9),
       pal = c("#ffffff", "#c7e9c0", "#74c476", "#238b45", "#00441b")),
  list(kind = "clifford", p = c(-1.8, -2.0, -0.5, -0.9),
       pal = c("#ffffff", "#dadaeb", "#9e9ac8", "#6a51a3", "#3f007d"))
)

ncol <- 2
nrow <- 2
cell_px <- 300

s <- vl_scene(6, 6, dpi = 100, bg = "white") |>
  push(vl_viewport(layout = grid_layout(
    widths  = vl_unit(rep(1, ncol), "null"),
    heights = vl_unit(rep(1, nrow), "null")
  )))

for (i in seq_along(gallery)) {
  a <- gallery[[i]]
  orbit <- attractor(a$kind, a$p)
  w <- window(orbit$x, orbit$y)
  img <- datashade(
    orbit$x, orbit$y,
    width = cell_px, height = cell_px,
    xlim = w$xlim, ylim = w$ylim,
    colors = a$pal, how = "eq_hist"
  )
  row <- (i - 1) %/% ncol + 1
  col <- (i - 1) %% ncol + 1
  s <- s |> push(vl_viewport(row = row, col = col)) |> draw(img) |> pop()
}

s
```

The shipped example (`inst/examples/attractors.R`) goes further: it draws twelve
panels with *random* parameters, keeping only orbits that fill a fair fraction
of the canvas (rejection sampling, since most random parameters collapse to a
point or diverge), each under a random colormap. Run it for a fresh gallery
every time:

```r
Rscript inst/examples/attractors.R gallery.png 1e7
```

## Recap

- `datashade()` bins a point cloud into a grid and colours cells by density, so
  cost tracks the grid size, not the point count, and overplotting becomes signal
  rather than noise.
- `how = "eq_hist"` keeps structure visible across orders of magnitude; `"log"`,
  `"cbrt"`, and `"linear"` are simpler alternatives.
- The result is one `raster_grob()`. Draw it inside a viewport whose scales match
  `xlim` / `ylim` to align it with axes, or tile many of them in a `grid_layout()`
  for a gallery.

## In a grammar

`datashade()` is the low-level engine: you supply the coordinates, the limits,
and the ramp. A grammar layer on top can wire all of that up from a plot spec.
[vellumplot](https://r-vellum.github.io/vellumplot/) exposes it as `mark_datashade()`,
which bins straight from a data frame and fits into a normal plot with scales,
guides, and facets. Its
[datashading article](https://r-vellum.github.io/vellumplot/articles/datashading.html)
pushes the same engine to its limit, shading the full US Census (about 306
million points) two ways.
```