Skip to content

Instantly share code, notes, and snippets.

@steveharoz
Last active August 29, 2015 14:02
Show Gist options
  • Save steveharoz/172ed1b1825f101be0f8 to your computer and use it in GitHub Desktop.
Save steveharoz/172ed1b1825f101be0f8 to your computer and use it in GitHub Desktop.
Mandelbrot image in R
library(ggplot2)
library(grid)
library(doParallel)
library(itertools)
resolution = 1000
iterations = 100
threshold = 10
# image range
data = expand.grid(
a = seq(-.65, -.5, length.out=resolution),
b = seq(.55, .70, length.out=resolution))
# initialize
data$ab = complex(real=data$a, imaginary=data$b)
data$d = iterations
# functions
distance = function(x, y) sqrt(x^2 + y^2)
mandelbrot = function(z, c) z^2 + c
# compute
compute = function (ab) {
z = ab
d = iterations # number of iterations until the point reaches a threashold distance
for (i in 1:iterations) {
z = mandelbrot(z, ab)
real = Re(z)
imaginary = Im(z)
dist = distance(real, imaginary)
d = ifelse(dist>threshold & d>i, i, d)
}
return(d)
}
# chunked parallel compute
computeChunk = function(x) {
# Create a cluster with 1 fewer cores than are available. Adjust as necessary
cl = makeCluster(detectCores() - 1) # 8-1 cores
registerDoParallel(cl)
# foreach
out = foreach(xc=isplitVector(x, chunks=getDoParWorkers()),
.export=c('mandelbrot', 'distance', 'iterations', 'threshold', 'compute'),
.combine='c') %dopar% {
compute(xc)
}
# finish
stopCluster(cl)
return(out)
}
# single threaded
#system.time(data$d <- compute(data$ab))
# multi-threaded
system.time(data$d <- computeChunk(data$ab))
# color gradient
colors = c('#000000', '#000000', '#cccc22', '#000000', '#ff0000', '#000000', '#112266', '#aaaa22', '#ffffff')
values = c(0, 0.1, 0.49, .491, .51, .52, .53, .99, 1)
plot = ggplot(data) +
geom_raster(aes(x=a, y=b, fill=log(d))) +
scale_fill_gradientn(colours=colors, values=values, guide='none') +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0))
# strip out everything but the plot
gt <- ggplot_gtable(ggplot_build(plot))
ge <- subset(gt$layout, name == "panel")
grid.draw(gt[ge$t:ge$b, ge$l:ge$r])
# save image
png('mandelbrot.png', width=resolution, height=resolution)
grid.draw(gt[ge$t:ge$b, ge$l:ge$r])
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment