mitchellgritts.com

The single most important quality in a piece of software is simplicity. It’s more important than doing the task you set out to achieve. It’s more important than performance. The reason is straightforward: if your solution is not simple, it will not be correct or fast. Drew DeVault

I've fallen victim to my own cleverness more times than I would like to admit. After reading Drew's post I've realized how often this happens to me. Let me show you an example.

Setup

The code below will get us setup to run the examples. I've also included some sf code to confirm the area of the polygon is indeed 1.

# calculate the area of a polygon, 
# in this case a square with an area of 1
x <- c(0, 1, 1, 0, 0)
y <- c(0, 0, 1, 1, 0)

# use sf to confirm area is indeed 1
sf::st_area(sf::st_polygon(list(cbind(x, y))))

# create vectors to traverse the sides of the polygon in
# r using the shoelace formula
x0 <- x[1:4]
x1 <- x[2:5]
y0 <- y[1:4]
y1 <- y[2:5]

The last few lines will be used to traverse the vertices of the polygon using the shoelace formula.

The first attempt

My first attempt at implementing the shoelace formula went as follows.

(Reduce(`+`, (x0 * y1)) - Reduce(`+`, (x1 * y0))) / 2
#> 1

The code does return the correct answer. I did account for vectorized operations in R. In other languages we would need to loop over the vector (x0, y1, ect...) in order to multiply them together. But is it simple? No it isn't.

I probably ended up at this implementation because I had just finished writing a lot of JavaScript. I tend to use reduce often in JavaScript, and in R Reduce is similar. I ended up with code that was half JavaScript inspired and half R inspired.

Encapsulation

Before moving on let's write a function to encapsulate our area function.

get_area <- function (polygon_matrix) {
  n <- nrow(polygon_matrix)
  
  x0 <- polygon_matrix[1:n - 1, 1]
  x1 <- polygon_matrix[2:n, 1]
  y0 <- polygon_matrix[1:n - 1, 2]
  y1 <- polygon_matrix[2:n, 2]

  (Reduce(`+`, (x0 * y1)) - Reduce(`+`, (x1 * y0))) / 2
}

Moving forward changes will be made to this function. And to be a little more rigorous, let's add some testing using the assertthat library. This will compare the results of sf::st_area to our get_area function. If assertthat::are_equal returns FALSE, our function isn't working properly.

polygon_matrix <- cbind(x, y)

# test
assertthat::are_equal(
  get_area(polygon_matrix),
  sf::st_area(sf::st_polygon(list(polygon_matrix)))
)

Attempt 2

Debugging is twice as hard as writing the code in the first place. Therefore, if you write the code as cleverly as possible, you are, by definition, not smart enough to debug it. Brian Kernighan

I was trying to be clever. I had found a perfect use case for the Reduce function. And not just once, but twice, in the same line of code! In fact, I was pleased with the cleverness of my code.

The hubris was short lived. In my cleverness I had overlooked the existance of the function sum. I couldn't believe how stupid I had been (my inner dialogue can be harsh).

So I rewrote the code using sum.

get_area <- function (polygon_matrix) {
  n <- nrow(polygon_matrix)
  
  x0 <- polygon_matrix[1:n - 1, 1]
  x1 <- polygon_matrix[2:n, 1]
  y0 <- polygon_matrix[1:n - 1, 2]
  y1 <- polygon_matrix[2:n, 2]

  (sum(x0 * y1) - sum(x1 * y0)) / 2
}

assertthat::are_equal(
  get_area(polygon_matrix),
  sf::st_area(sf::st_polygon(list(polygon_matrix)))
)
#> [1] TRUE

Indeed, it returns the correct answer. The main logic in the function is still one line of code. It is simple. But it can be simplified even more!

Attempt 3

Alas, it can be even simpler. And potentially faster. sum is used twice. sum is vectorized and I've made use of it in the function above. But I forgot about the order of operations, and the actual mathematical notation of the shoelace formula.

Back to the drawing board.

get_area <- function (polygon_matrix) {
  n <- nrow(polygon_matrix)
  
  x0 <- polygon_matrix[1:n - 1, 1]
  x1 <- polygon_matrix[2:n, 1]
  y0 <- polygon_matrix[1:n - 1, 2]
  y1 <- polygon_matrix[2:n, 2]

  sum(x0 * y1 - x1 * y0) / 2
}

assertthat::are_equal(
  get_area(polygon_matrix),
  sf::st_area(sf::st_polygon(list(polygon_matrix)))
)

Finally, correct and simple. And since there is only a single sum function call it should be faster. I've not performed any benchmarking to confirm this.

Summary

Just to be certain, let's try this code on a more complex polygon to make sure it works. The area of the polygon is 30. Let's see if the get_area function works properly.

x <- c(3, 5, 9, 12, 5, 3)
y <- c(4, 6, 5, 8, 11, 4)

polygon_matrix <- cbind(x, y)

assertthat::are_equal(
  get_area(polygon_matrix),
  sf::st_area(sf::st_polygon(list(polygon_matrix)))
)
#> [1] TRUE

The first iteration of the get_area function was complex. Two needless calls to a function that isn't very common (Reduce), using a syntax that isn't super familiar (using an operator within the back ticks).

The second iteration doesn't use Reduce, but two redundant instances of the sum function. Simpler, but not the simplest. The final implementation is indeed the simplest, and it's proven to be correct (and I think faster).