Category Archives: Fractals

The Sound Of Mandelbrot Set

Music is the pleasure the human soul experiences from counting without being aware that it is counting (Gottfried Leibniz)

I like the concept of sonification: translating data into sounds. There is a huge amount of contents in the Internet about this technique and there are several packages in R to help you to sonificate your data. Maybe one of the most accessible is tuneR, the one I choosed for this experiment. Do not forget to have a look to playitbyr: a package that allows you to listen to a data.frame in R by mapping columns onto sonic parameters, creating an auditory graph, as you can find in its website. It has a very similar syntaxis to ggplot. I will try to post something about playitbyr in the future.

Let me start plotting the Mandelbrot Set. I know you have seen it lot of times but it is very easy to plot in with R and result is extremely beautiful. Here you have four images corresponding to 12, 13, 14 and 15 iterations of the set’s generator. I like a lot how the dark blue halo around the Set evaporates as number of iterations increases.

Mandelbrot12Mandelbrot13Mandelbrot14Mandelbrot15

And here you have the Set generated by 50 iterations. This is the main ingredient of the experiment:Mandelbrot50

Mandelbrot Set is generated by the recursive formula xt+1=xt2+c, with x0=0. A complex number c belongs to the Mandelbrot Set if its module after infinite iterations is finite. It is not possible to iterate a infinite number of times so every representation of Mandelbrot Set is just an approximation for a usually big amount of iterations. First image of Mandelbrot Set was generated in 1978 by Robert W. Brooks and Peter Matelski. You can find it here. I do not know how long it took to obtain it but you will spend only a couple of minutes to generate the ones you have seen before. It is amazing how computers have changed in this time!

This iterative equation is diabolical. To see just how pathological is, I transformed the succession of modules of xt generated by a given c in a succession of sounds. Since it is known that if one of this iterated complex numbers exceeds 2 in module then it is not in the Mandelbrot Set, frequencies of these sounds are bounded between 280 Hz (when module is equal to zero) and 1046 Hz (when module is equal or greater to 2). I called this function CreateSound. Besides the initial complex, you can choose how many notes and how long you want for your composition.

I tried with lot of numbers and results are funny. I want to stand out three examples from the rest:

  • -1+0i gives the sequence 0, −1, 0, −1, 0 … which is bounded. Translated into music it sounds like an ambulance siren.
  • -0.1528+1.0397i that is one of the generalized Feigenbaum points, around the Mandelbrot Set is conjetured to be self-similar. It sounds as a kind of Greek tonoi.
  • -3/4+0.01i which presents a crazy slow divergence. I wrote a post some weeks ago about this special numbers around the neck of Mandelbrot Set and its relationship with PI.

All examples are ten seconds length. Take care with the size of the WAV file when you increase duration. You can create your own music files with the code below. If you want to download my example files, you can do it here. If you discover something interesting, please let me know.

Enjoy the music of Mandelbrot:

# Load Libraries
library(ggplot2)
library(reshape)
library(tuneR)
rm(list=ls())
# Create a grid of complex numbers
c.points <- outer(seq(-2.5, 1, by = 0.002),1i*seq(-1.5, 1.5, by = 0.002),'+')
z <- 0
for (k in 1:50) z <- z^2+c.points # Iterations of fractal's formula
c.points <- data.frame(melt(c.points))
colnames(c.points) <- c("r.id", "c.id", "point")
z.points <- data.frame(melt(z))
colnames(z.points) <- c("r.id", "c.id", "z.point")
mandelbrot <- merge(c.points, z.points, by=c("r.id","c.id")) # Mandelbrot Set
# Plotting only finite-module numbers
ggplot(mandelbrot[is.finite(-abs(mandelbrot$z.point)), ], aes(Re(point), Im(point), fill=exp(-abs(z.point))))+
geom_tile()+theme(legend.position="none", axis.title.x = element_blank(), axis.title.y = element_blank())
#####################################################################################
# Function to translate numbers (complex modules) into sounds between 2 frequencies
#   the higher the module is, the lower the frequencie is
#   modules greater than 2 all have same frequencie equal to low.freq
#   module equal to 0 have high.freq
#####################################################################################
Module2Sound <- function (x, low.freq, high.freq)
  {
    if(x>2 | is.nan(x)) {low.freq} else {x*(low.freq-high.freq)/2+high.freq}
  } 
#####################################################################################
# Function to create wave. Parameters:
#    complex     : complex number to test
#    number.notes: number of notes to create (notes = iterations)
#    tot.duration.secs: Duration of the wave in seconds
#####################################################################################
CreateSound <- function(complex, number.notes, tot.duration.secs)
{
  dur <- tot.duration.secs/number.notes
  sep1 <- paste(", bit = 16, duration= ",dur, ", xunit = 'time'),sine(")
  sep2 <- paste(", bit = 16, duration =",dur,",  xunit = 'time'))")
  v.sounds <- c()
  z <- 0
  for (k in 1:number.notes) 
  {
    z <- z^2+complex
    v.sounds <- c(v.sounds, abs(z))
  }
  v.freqs <- as.vector(apply(data.frame(v.sounds), 1, FUN=Module2Sound, low.freq=280, high.freq=1046))
  eval(parse(text=paste("bind(sine(", paste(as.vector(v.freqs), collapse = sep1), sep2)))
}
sound1 <- CreateSound(-3/4+0.01i     , 400 , 10) # Slow Divergence
sound2 <- CreateSound(-0.1528+1.0397i, 30  , 10) # Feigenbaum Point
sound3 <- CreateSound(-1+0i          , 20  , 10) # Ambulance Siren
writeWave(sound1, 'SlowDivergence.wav')
writeWave(sound2, 'FeigenbaumPoint.wav')
writeWave(sound3, 'AmbulanceSiren.wav')

I Need A New Computer To Draw Fractals!

Computer Science is no more about computers than astronomy is about telescopes (E. W. Dijkstra)

Some days ago I published a post about how to build fractals with R using Multiple Reduction Copt Machine (MRCM) algorithm. Is that case I used a feature of the grid package that allows you to locate objects easily into the viewPort avoiding to work with coordinates. It does not work well if you want to divide your seed image into five subimages located in the vertex of a regular pentagon. No problem: after refreshing some trigonometric formulas and after understanding how to work with coordinates I felt strong enough to program the Final-MRCM-Fractal-Builder. But here comes the harsh reality. My computer crashes when I try to go beyond five degrees of depth. Imposible. In the example of Sierpinski’s triangle, where every square in divided into three small ones, I reached seven degrees of depth. I am deeply frustrated. These are drawings for 1, 2, 3 and 5 degrees of depth.

Rplot01c Rplot02c Rplot03cRplot05c

Please, if someone modifies code to make it more efficient, let me know. I used circles in this case instead squares. Here you have it:

library(grid)
grid.newpage()
rm(list = ls())
ratio <- 0.4
pmax <- 5 # Depth
vp1 <- viewport(w=1, h=1)
vp2 <- viewport(w=ratio, h=ratio, just=c(0.75*sin(2*pi*1/5)+0.5, 0.75*cos(2*pi*1/5)+0.75*pi*1/5))
vp3 <- viewport(w=ratio, h=ratio, just=c(0.75*sin(2*pi*0/5)+0.5, 0.75*cos(2*pi*0/5)+0.75*pi*1/5))
vp4 <- viewport(w=ratio, h=ratio, just=c(0.75*sin(2*pi*2/5)+0.5, 0.75*cos(2*pi*2/5)+0.75*pi*1/5))
vp5 <- viewport(w=ratio, h=ratio, just=c(0.75*sin(2*pi*3/5)+0.5, 0.75*cos(2*pi*3/5)+0.75*pi*1/5))
vp6 <- viewport(w=ratio, h=ratio, just=c(0.75*sin(2*pi*4/5)+0.5, 0.75*cos(2*pi*4/5)+0.75*pi*1/5))
pushViewport(vp1)
grid.rect(gp=gpar(fill="white", col=NA))
m <- as.matrix(expand.grid(rep(list(2:6), pmax)))
for (j in 1:nrow(m))
{
for(k in 1:ncol(m)) {pushViewport(get(paste("vp",m[j,k],sep="")))}
grid.circle(gp=gpar(col="dark grey", lty="solid", fill=rgb(sample(0:255, 1),sample(0:255, 1),sample(0:255, 1), alpha= 95, max=255)))
upViewport(pmax)
}

What The Hell Is Pi Doing Here?

Nothing in Nature is random … A thing appears random only through the incompleteness of our knowledge (Benedict Spinoza)

This is one of my favorite mathematical mysteries. In 1991 David Boll was trying to confirm that the neck of the Mandelbrot Set is 0 in thickness. Neck is located at -0.75+0i (where two biggest circles meet each other). Mandelbrot50He tried with complex numbers like -0.75+εi for small values of ε demonstrating the divergence of all these numbers. And here comes the mystery: multiplying ε and the corresponding number of iterations it took for the iterate to diverge, gives an approximation of π that is within ±ε. Is not fascinating? I replicated David Boll’s experiment for positive and negative values of ε. I draw results as follows:

mandelbrot

Before doing it, I thought I was going to find some pattern in the graphic. Apart from the mirror effect produced by the sign of ε, there is nothing recognizable. Convergence is chaotic. Here you have the code. This example is also nice to practice with ggplot2 package, one of the totems of R:

i<-0    # Counter of iterations
x  <- 0 # Initialization
while (Mod(x) <= 2)
{
x <- x^2+(c+complex(real = 0, imaginary = e))
i <- i+1
}
i
}
results <- as.data.frame(c(NULL,NULL))
for (j in 1:length(epsilons))
{results <- rbind(results, c(epsilons[j], testMSConvergence(epsilons[j])))}
colnames(results) <- c('epsilon', 'iterations')
dev.off()
p <- ggplot(results, aes(epsilon,abs(epsilon)*iterations))+
xlab("epsilon")+
ylab("abs(epsilon)*iterations")+
opts(axis.title.x=theme_text(size=16)) +
opts(axis.title.y=theme_text(size=16)) +
ggtitle("How to Estimate Pi Using Mandelbrot Set's Neck")+
theme(plot.title = element_text(size=20, face="bold"))
p <- p + geom_ribbon(data=results,aes(ymin=abs(epsilon)*iterations-abs(epsilon),ymax=abs(epsilon)*iterations+abs(epsilon)), alpha=0.3)
p <- p + geom_abline(intercept = pi, , slope = 0, size = 0.4, linetype=2, colour = "black", alpha=0.8)
p <- p + geom_line(colour = "dark blue", size = 1, linetype = 1)
p <- p + geom_text(x = 0, y = pi, label="pi", vjust=2, colour="dark blue")
p <- p + geom_point(x = 0, y = pi, size = 6, colour="dark blue")
p + geom_point(x = 0, y = pi, size = 4, colour="white")

Building Affine Transformation Fractals With R

Clouds are not spheres, mountains are not cones, coastlines are not circles and bark is not smooth, nor does lightning travel in a straight line (Benoit Maldelbrot)

Fractals are beautiful, hypnotics, mysterious. Cantor set has as many points as the real number line but has zero measure. After 100 steps, the Koch curve created from a 1 inch segment is long enough  to wrap around the Earth at the equator nearly four thousand times. The Peano Curve is a line that has the same dimension as a plane. Fractals are weird mathematical objects. Fractals are very cool.

One way to build fractals is using the Multiple Reduction Copy Machine (MRCM) algorithm which uses affine linear transformations over a seed image to build fractals. MRCM are iterative algorithms that perform some sort of copy+paste task. The idea is quite simple: take a seed image, transform it (clonation, scalation, rotation), obtain the new image and iterate.

To create the Sierpinsky Gasket Fractal you part from a square. Then you divide it into three smaller squares, locate them as a pyramid and iterate doing the same with avery new square created. Making these things is very easy with grid package. Defining the division (i.e. the affine transformation) properly using viewPort function and navigating between them is all you need. Here you have the Sierpinsky Gasket Fractal with 1, 3, 5 and 7 levels of depth. I filled in squares with random colours (I like giving some touch of randomness to pictures). Here you have pictures:

Rplot01 Rplot03

Rplot05 Rplot07

And here you have the code. Feel free to build your own fractals.

library(grid)
rm(list = ls())
grid.newpage()
pmax <- 5 # Depth of the fractal
vp1=viewport(x=0.5,y=0.5,w=1, h=1)
vp2=viewport(w=0.5, h=0.5, just=c("centre", "bottom"))
vp3=viewport(w=0.5, h=0.5, just=c("left", "top"))
vp4=viewport(w=0.5, h=0.5, just=c("right", "top"))
pushViewport(vp1)
m <- as.matrix(expand.grid(rep(list(2:4), pmax)))
for (j in 1:nrow(m))
{
for(k in 1:ncol(m)) {pushViewport(get(paste("vp",m[j,k],sep="")))}
grid.rect(gp=gpar(col="dark grey", lty="solid",
fill=rgb(sample(0:255, 1),sample(0:255, 1),sample(0:255, 1), alpha= 95, max=255)))
upViewport(pmax)
}