Una lección de geometría

Necesitaremos datos para realizar el benchmark, y ya que hablamos de polígonos, sería conveniente contar con una función que los genere aleatoriamente. Recordemos que un polígono es una figura plana finita, delimitada por una cadena cerrada (el límite) de segmentos lineales (los lados). Si el límite no se interseca a sí mismo, entonces el polígono es simple. Para generalizar la nomenclatura nos referiremos al polígono de n lados como un “n-gon”. Ahora, aunque no lo parezca, generar polígonos completamente aleatorios es un problema más o menos complicado. El siguiente método simplifica considerablemente este problema:

1. Ubicar n puntos aleatorios sobre una circunferencia de radio definido.
2. Ordenar los puntos en sentido antihorario, comenzando por cualquiera de ellos.
3. Conectar los puntos en dicho orden, cerrando y rellenando el polígono.

Debido a que todos sus vértices yacen sobre la circunferencia circunscrita, el polígono así generado será cíclico. Pero no todos los polígonos poseen una circunferencia circunscrita; los siguientes ejemplos -no obtenibles con este método- lo demuestran:

Hasta aquí, los polígonos presentados, cíclicos o no, son todos convexos; significa que todos sus ángulos interiores son menores a 180º. Todos los polígonos cíclicos son convexos, pero lo contrario no es verdad. La principal propiedad de la convexidad es que, para dos puntos cualesquiera del polígono (en su interior o en su límite), el segmento que los conecta nunca sale al exterior del polígono. En contraste, los polígonos cóncavos (aquellos con al menos un ángulo interior mayor a 180º), carecen de esta propiedad. Una manera sencilla (aunque no infalible) de obtener polígonos cóncavos es insertar este paso en el método anterior:

2.5. Mover aleatoriamente los puntos hacia el centro de la circunferencia, con una cantidad no mayor al radio de la misma.

Dependiendo de la magnitud del acercamiento de los vértices al centro, los polígonos podrían o no dejar de ser convexos, pero siempre dejarán de ser cíclicos. Además, siempre serán estrellados. Un polígono estrellado posee al menos un punto en su interior desde el cual es posible observar todo el límite; en este método, al menos el centro de la circunferencia cumplirá esta propiedad. Todos los polígonos convexos son estrellados, pero no todos los cóncavos lo son; los siguientes ejemplos son cóncavos, pero no estrellados:

En definitiva, este método es una simplificación porque solamente puede producir polígonos estrellados; sin embargo, para nuestro propósito, es suficientemente simple y aleatorio. Con esto concluye la lección de geometría.

Un generador de polígonos

La siguiente función ngon utiliza el método ya descrito para generar aleatoriamente un polígono estrellado de n lados, convexo o cóncavo, contenido en una circunferencia de radius específico centrada en el origen.

ngon = function(n, min_radius = NA, radius = 1){

   stopifnot(is.integer(n), n > 2L)
   stopifnot(is.numeric(radius), radius > 0)
   
   if(!is.na(min_radius)){
      stopifnot(is.double(min_radius), min_radius > 0, min_radius < 1)
      radius = radius * runif(n, min_radius, 1)
   }

   theta = sort(runif(n, max = 2 * pi))
   ngon = radius * matrix(c(cos(theta), sin(theta)), n)
   colnames(ngon) = c("x", "y")
   return(`class<-`(ngon, c("ngon", "matrix")))
}

El resultado predeterminado de ngon sería un polígono cíclico, con circunferencia circunscrita de radio unitario. La probabilidad de obtener un polígono cóncavo deriva del argumento min_radius que permite acercar los vértices al centro, hasta cierto radio mínimo (como porcentaje del radio de la circunferencia). Tanto este acercamiento, como la ubicación de los vértices en la circunferencia, son calculados con una distribución uniforme, vía runif.

Para facilitar el trabajo más adelante, se han definido dos métodos¹ para los objetos de clase ngon (uno para graficar y otro para transformar a la clase units), y una función para transformar a la clase sfg.

plot.ngon = function(x, pal = tmaptools::get_brewer_pal("Dark2", 15, plot=FALSE), ...){
   plot.default(x, asp = 1, bty = "n", type = "n", xaxt = "n", yaxt = "n", ...)
   polygon(x, col = sample(pal, 1), border = NA)
}

library(units)

set_units.ngon = function(x, value) units::set_units(
   `class<-`(x, "matrix"), value, mode = "standard")

library(sf)

as_polygon = function(x) st_polygon(list(x[c(1:nrow(x), 1), ]))

set.seed(57); par(mfrow = c(2, 3)) # la "semilla" determina cuáles polígonos obtendré

for(i in c(3:4L, 6L, 10L, 30L, 90L)) plot(ngon(i), xlab = paste0(i, "-gon"), ylab = "")

ngon(1e4L) %>% set_units("cm") %>% shoelace_area() %>% format(digits = 8)

## [1] "3.1415915 [cm^2]"

set.seed(32); par(mfrow = c(3, 3))

for(i in c(5L, 10L, 30L)){
   for(j in c(.7, .4, .1)){
      plot(ngon(i, j), xlab = paste0(i, "-gon"), ylab = paste("min", j))
   }
}

Efectuar benchmark

Finalmente poseemos los datos y el conocimiento necesarios para evaluar la velocidad de ejecución de dos funciones que calculan áreas. En pocas palabras, el benchmark consiste en ejecutar cierta expresión (una función o algún fragmento pequeño de código) muchas veces, midiendo cada vez el tiempo invertido en ello. Obviamente el tiempo dependerá de la capacidad del procesador utilizado; sin embargo, la utilidad del benchmark reside en ejecutar no una, sino varias expresiones similares. Comparando los tiempos de ejecución se obtiene un criterio -que algunas veces resulta sufiente- para decidir cuál expresión es la mejor resolviendo cierto problema.

De regreso a la definición de ngon, existen tres parámetros con los que se puede jugar y evaluar si ocurre algún cambio. Sin embargo, como se verá, el radio de la circunferencia circunscrita y la convexidad del polígono generado son irrelevantes; por esta razón se estudió su comportamiento con la función ngon solamente. Todos los benchmarks fueron ejecutados con microbenchmark.

library(microbenchmark)

for(i in c(10, 30, 100)) assign(paste0("r", i), ngon(100L, radius = 10 ^ i))

radBmark = microbenchmark(shoelace_area(r10), shoelace_area(r30), shoelace_area(r100))

radBmark = dplyr::filter(radBmark, time < 15e3)

## Unit: microseconds
##                 expr min  lq  mean median   uq  max neval
##   shoelace_area(r10) 8.2 8.5 9.020    8.9 9.30 11.8    99
##   shoelace_area(r30) 8.0 8.4 8.984    8.8 9.35 12.1    99
##  shoelace_area(r100) 8.2 8.5 9.027    8.9 9.30 11.3   100

for(i in c(9, 5, 1)) assign(paste0("r.", i), ngon(1000L, i / 10))

conBmark = microbenchmark(shoelace_area(r.9), shoelace_area(r.5), shoelace_area(r.1))

verBmark = setNames( , c(1e2L, 3e2L, 1e3L, 3e3L, 1e4L, 3e4L)) %>%
   lapply(function(ver) microbenchmark(
      shoelace_area(gon), st_area(sfg), times = 1e3,
      setup = { gon = ngon(ver); sfg = as_polygon(gon) }
      ))

Conclusiones

En este artículo se habló de algunos conceptos de geometría y de cómo implementarlos en código R, aparentemente con el propósito de generar polígonos y comparar dos funciones que permiten calcular sus áreas. Pero la verdad es que el propósito siempre fue aprender a realizar benchmarks. Desde el principio era previsible que sf::st_area tuviese una ejecución más larga, no solo porque debe verificar que los datos sean adecuados (de clase sf / sfc / sfg); también porque está preparada para tratar con los distintos tipos de geometrías existentes en el análisis espacial. En el siguiente ejemplo, el área de la geometría multipolygon no podría ser calculada con nuestra función.

data("NLD_prov", package = "tmap")

(NLD_prov$area = set_units(st_area(NLD_prov), "km2"))

## Units: [km^2]
##  [1] 2392 3518 2680 3405 1465 5118 1443 2857 3103 1834 5053 2210

NLD_cent = unlist(st_centroid(NLD_prov)$geometry)

plot(NLD_prov["area"], main = "Netherlands provinces' area [km²]", reset = FALSE)

text(NLD_cent[2*1:12-1], NLD_cent[2*1:12], sub("-", "\n", NLD_prov$name), col = "grey")

Adicionalmente, ambas fórmulas se ejecutan con gran velocidad: cuando se calculó el área de un 30000-gon, el tiempo promedio fue algo mayor a mil microsegundos (eso es la milésima parte de un segundo). Una duda final: ¿en sf::st_area se aplica la fórmula shoelace, o existe un mejor algoritmo? Observando los tiempos de nuestra shoelace_area, no obstante, parece que sí es la solución más eficiente.