Happy Pi Day everyone!!

Overview: Today is 3-14, and those digits are the first three digits of Pi. To celebrate, I estimated Pi using a lightsaber battle in R. This is a fun math demonstration where we first build an army of saber-wielding critters in the park and then cut transects across the battlefield and count the number of contacts we make with lightsabers along the way. After the battle, we take a weighted ratio of hits vs. misses to estimate Pi.

Our adventure begins in a futuristic world where urban wildlife have aquired lightsabers and taken over Montreal.

Figure 1: Future squirels battle for food scraps near homeless man’s dome tent.

We are about to join the battle. It is probably a good time to learn how lightsabers work. We represent a lightsaber with a line, meaning that we assume the saber has no width. We describe that line using xy coordinates in a matrix.

set.seed(1.1)
light_saber_length <- 20
number_of_animal_sabers_in_a_park <- 5000
distance_between_transects <- 40.1

point_1 <- c(-(light_saber_length/2),(light_saber_length/2))
point_2 <- c(0,0)
single_saber <-cbind(point_1, point_2)
rownames(single_saber)<- c("point_1","point_2")
colnames(single_saber)<- c("x","y")
single_saber

##           x y
## point_1 -10 0
## point_2  10 0

Figure 2: Represent our lightsaber as a line with no width. A single beam of light.

To move the object, we multiply that matrix by transformation matrices to rotate and move the object around. We rotate the object by multiplying the object matrix by the rotation matrix specifying a random angle.

  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_saber)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% t(matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2))
  
after_rotation

##              [,1]      [,2]
## point_1  1.100398 -9.939272
## point_2 -1.100398  9.939272

Figure 3: We rotate the lightsaber around the origin to change its angle.

We then move the saber to a location the battlefield. Note: order is important. We must always rotate around the origin first before moving into place on the landscape. It does not work the other way.

  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(
    c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  colnames(new_coord) <- c("x","y")
  rownames(new_coord) <- c("Point_1","Point_2")

Figure 4: After the angle is set, translate the lightsaber to a random location on the battlefield.

Create a function to rotate and distribute lightsabers to all the critters on the battlefield.

critter_saber <- function(){
  single_saber <-cbind(c(-5,5),c(0,0))
  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_saber)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2, byrow=TRUE )
  
  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(
    c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  colnames(new_coord) <- c("x","y")
  rownames(new_coord) <- c("Point_1","Point_2")
  
  return(new_coord)   
}

critter_saber()

##                x        y
## Point_1 21.93899 84.59096
## Point_2 30.33012 79.15139

Repeat it many times to create an army of lightsaber-wielding critters. We assume that all the critters are battling vigorously and indiscriminately with each other so they are randomly distributed across space and the angles of all the sabers are uniformly distributed.

for(i in 1:number_of_animal_sabers_in_a_park){
critter_saber()
}

An army of lightsaber-wielding critters in the middle of the battlefield. This is the position of each critter's lightsaber at the moment we pass them on our transect.

Figure 5: An army of lightsaber-wielding critters in the middle of the battlefield. This is the position of each critter’s lightsaber at the moment we pass them on our transect.

Figure 6: We better get in there! Things aren’t looking good.

surface_line <- function(line_height = 0){
  matrix(c(0,100,line_height,line_height), 2, 2, byrow=FALSE)
  }

##      [,1] [,2]
## [1,]    0   20
## [2,]  100   20

Figure 7: My transect through the army of critters.

one_saber <- critter_saber()
one_surface_line <- surface_line(30)
line.line.intersection(one_saber[1,], one_saber[2,], one_surface_line[1,], 
                       one_surface_line[2,], interior.only = TRUE)

## [1] 75.73794 30.00000

Figure 8: Build a mechanism to identify when a lightsaber has crossed my transect line.

Figure 9: Team prepares for transects

Prepare a list of all transects and all critter sabers so we can calculate contacts between them

list_of_matches <- array(NA,dim = c(2,2,number_of_animal_sabers_in_a_park))
list_of_surface_lines <- array(NA,dim = c(2,2,6))

for(i in 1:number_of_animal_sabers_in_a_park){
  list_of_matches[,,i] <- critter_saber()
}

line_placement <- seq(0,100, by=20)
for(i in 1:6){
  list_of_surface_lines[,,i] <- surface_line(line_placement[i])
}

The battle is on! We rush along our transects and dual with every lightsaber that crosses our path. We are so skilled that we hit every lightsaber to cross our transect but never receive a hit or hit anything off of our transect line.

Figure 10: The battle is on. We charge across the battlefield and hit every critter saber along our transect.

The battle is over! How did we do? We better count up all our hits and compare them to our misses.

crosses <- matrix(NA, 1, 10, byrow=TRUE)

for(i in 1:dim(list_of_matches)[3]){
for(j in 1:dim(list_of_surface_lines)[3]){
  crosses <- rbind(crosses, c(
    line.line.intersection(list_of_matches[1,,i], list_of_matches[2,,i], 
                           list_of_surface_lines[1,,j],     
    list_of_surface_lines[2,,j], interior.only = TRUE), 
    
    list_of_matches[1,,i], list_of_matches[2,,i], list_of_surface_lines[1,,j], 
    list_of_surface_lines[2,,j])
    )
  }
}

crosses <- crosses[which(is.na(crosses[,1]) == FALSE),]
colnames(crosses) <- c("x.cross", "y.cross", "x.p1", "y.p1", "x.p2", 
                    "y.p2", "x.surface.p1", "y.surface.p1", "x.surface.p2", 
                    "y.surface.p2")
head(crosses)

##       x.cross y.cross     x.p1     y.p1     x.p2     y.p2 x.surface.p1
## [1,] 63.08661      20 64.26028 15.53502 61.71804 25.20647            0
## [2,] 14.37017      80 12.61642 83.36074 17.24275 74.49524            0
## [3,] 49.44207      80 44.92515 85.23335 51.45904 77.66313            0
## [4,] 47.22221      20 49.37535 17.30391 43.13500 25.11787            0
## [5,] 24.43356      20 23.31048 15.67796 25.82545 25.35654            0
## [6,] 56.12789      20 53.77911 24.62028 58.31078 15.70602            0
##      y.surface.p1 x.surface.p2 y.surface.p2
## [1,]           20          100           20
## [2,]           80          100           80
## [3,]           80          100           80
## [4,]           20          100           20
## [5,]           20          100           20
## [6,]           20          100           20

Figure 11: All the lightsabers we struck on our transects.

We take the ratio of the lightsabers we made contact with against the ones we missed, scaled by distance, and we get a great approximation of Pi!

pi_estimate <- (2* light_saber_length * number_of_animal_sabers_in_a_park)/(distance_between_transects * length(crosses[,1]))
options(digits=20)
pi_estimate

## [1] 3.1427417593382642735

Estimating Pi using a lightsaber battle