The title of chapter 5 in my Guerrilla Capacity Planning book is, “Evaluating Scalability Parameters,” and underneath it you’ll see this quote:

“With four parameters I can fit an elephant. With five I can make his trunk wiggle.” —John von Neumann

In that vein, Guerrilla alumnus Stephen O’C. pointed me at a recent blog post and paper (PDF) that draws an elephantine curve using just 4 fitting coefficients or parameters. Stephen also sent me his translation of the Python code into R. Previous efforts apparently had required some 30 parameters. The secret to the success of this latest example is plotting the elephant in the complex plane by summing certain Fourier modes. That’s all very cool but I was surprised to see that the output was static (no wiggles), even though 5 parameters are defined. That shortcoming, however, provided me with the impetus to try out R’s animation package and here’s the result.

Notice that my elephant not only wiggles his trunk but he also winks—a wiggling winking pink elephant. Actually, I think he looks more like a winking woolly mammoth. 🙂

Since I also know a thing or two about complex numbers and Fourier transforms, I found the Python code a bit obtuse but I didn’t want to spend a lot of time rewriting it. Instead, I just looked for opportunities to:

• use the R animation package
  
• reorganize the code
  
• use R constructs such as lapply or sapply instead of loops
  
• use the fft (fast fourier transform) function in R
  

and that’s what I discuss in the remainder of this post. First, let’s consider a little bit about how this thing works.

The resulting shape (and wiggling) of this elephant is controlled by a set of 5 (complex) parameters and their rotated complement. That’s actually 20 numbers, but who’s counting?

require(MASS)
library(animation)

mkmovie = TRUE  #guarantees some form of output

param <- c(50-30i,18+8i,12-10i,-14-60i,1+20i)
parar <- param * exp(1i*pi/2)  #rotate by 90 degrees

If we apply each of the first 4 parameters in succession, the results look like this.

The first parameter produces an ellipse, the second makes a dent in the ellipse, while all 4 parameters together produce the base elephant shape. Roughly speaking, each parameter (p) produces p + 1 lobes in the 2D figure ... or are they limbs? The following R function assigns the parameters to the Fourier coefficients (Cx and Cy) to draw "pinky" the elephant.

pinky <- function() {
 Cx <- as.complex(rep(0,length(param)))
 Cy <- as.complex(rep(0,length(param))) 
 tv <- seq(0,2*pi,length=1000)

 for (i in 1:2) { #movie frames
  Cx[1] <- parar[1] + Im(param[1])
  Cx[2] <- parar[2] + Im(param[2])
  Cx[3] <- Re(param[3])
  Cx[4] <- Re(param[5]) - (i-1)
  Cx[5] <- Re(param[4])
  
  Cy[1] <- param[1] - Re(param[1]) + Im(param[4])
  Cy[2] <- param[2] - Re(param[2])
  Cy[3] <- param[3] - Re(param[3])
  
  x <- c(fourier(tv, Cx))
  y <- c(fourier(tv, Cy))
  
  plot(y, -x, type="l", col='red', lwd=10, axes=FALSE, ylab='', xlab='')
  lines(y, -x, type="l", col='pink', lwd=4)
  if (i > 1) points(Im(param[5]), Im(param[5]), col='black', pch=126, cex=2)
     else points(Im(param[5]), Im(param[5]), col='black', pch=20, cex=2)
  }
}

Other combinations of these parameters can produce such things as Lissajous figures, just like I used to make on my father's oscilloscope (when he wasn't using it to fix our TV). I also grew up seeing another Lissajous figure: the logo of the Australian Broadcasting Corporation.

The vectors Cx and Cy are used in the Fourier summation, which I've written as a recursive function.

fourier <- function(tt,cc) {
 wt <- rep(0, length(tt)) 
 fsum <- function(n) {
   if (n > 0) wt <- wt + fsum(n-1) + Re(cc[n]) * cos(n*tt) + Im(cc[n]) * sin(n*tt)
   return(wt)
 }
 fsum(length(cc))
}

This tells me immediately that there's no point using fft in R because there, the purpose is to determine the parameter values (as outputs) that produce a given periodic function. Here, the parameter values are being chosen deliberately as inputs to draw a function that resembles an elephant. That's also antithetical to von Neumann's statement but let's not get too pedentic.

I also tried to find a way to use lapply directly but I ran into a conflict between the size of the tt (time) vector versus the number of summation terms. In any event, I was quickly facing the prospect of writing some additional functions for lapply and it all started to feel like I was overloading its intent. Worse yet, I would also have ended up with more code than the original for-loop. Ultimately, I just replaced it with the recursive routine above and declared victory.

Finally, the two frames from the pinky function can be output optionally as an animated GIF, such as I've used here, or as two static plots in the standard R graphics window.

if (mkmovie) {
 aopt = ani.options(interval = 0, nmax = 2)
 saveMovie(pinky(), interval = 0.25, outdir = getwd(), width = 400, height = 400)
 ani.options(aopt)
} else pinky()

Who knew that elephants could be so mathematically interesting. Animation is also a powerful tool for performance data visualization and I'll say more about that in the upcoming GDAT class.

Thanks, Stephen!