マンデルブロ集合をRで計算する。並列化の練習もしてみる。
基本的にはMy Christmas Gift: Mandelbrot Set Computation In Python、How To Quickly Compute The Mandelbrot Set In Pythonをなぞっている。GPUによる並列化はまだ出来ていない。
定義
wikipediaより
漸化式
で定義される複素数列$z_{n}\in \mathbb{C}$が$n\to\infty$の極限で無限大に発散しない複素数$c$全体が作る集合を、マンデルブロ集合という(以下$\mathbb{M}$と表す)。
性質
単純なRコード
mandelbrot_1 <- function(c, max.iter=128) { z <- c for (rep in 1:max.iter) { if (Mod(z) > 2) return(list(z=z, k=rep-1)) z <- z*z + c } return(list(z=z, k=max.iter)) } mandelbrot_naive <- function(xmin=-2.0, xmax=0.5, ymin=-1.25, ymax=1.25, nx=500, ny=500, max.iter=128) { X <- seq(xmin, xmax, length.out=nx) Y <- seq(ymin, ymax, length.out=ny) Z <- matrix(0.0, nrow=nx, ncol=ny) K <- matrix(0.0, nrow=nx, ncol=ny) for (i in 1:nx) for (j in 1:ny) { tmp <- mandelbrot_1(X[i] + 1i*Y[j], max.iter) Z[i, j] <- tmp$z K[i, j] <- tmp$k } return(list(X=X, Y=Y, Z=Z, K=K)) } xmin=-2.0; xmax=0.5; ymin=-1.25; ymax=1.25; nx=500; ny=500; max.iter=128 M <- mandelbrot_naive(xmin, xmax, ymin, ymax, nx, ny, max.iter) mandelbrot_plot <- function(M) { cols <- c( colorRampPalette(c("#00008b", "blue", "#87cefa", "#87ceeb", "#b0e0e6", "white"))(100), "black") image(M$X, M$Y, M$K, col=cols, xlab="Re(c)", ylab="Im(c)") } mandelbrot_plot(M)
ベクトル化したRコード(1)
mandelbrot_vectorized_1 <- function(xmin=-2.0, xmax=0.5, ymin=-1.25, ymax=1.25, nx=500, ny=500, max.iter=128) { X <- seq(xmin, xmax, length.out=nx) Y <- seq(ymin, ymax, length.out=ny) C <- outer(X, 1i*Y, FUN="+") Z <- C K <- matrix(0.0, nrow=nx, ncol=ny) for (rep in 1:max.iter) { index <- which(Mod(Z) <= 2) Z[index] <- Z[index]^2 + C[index] K[index] <- K[index] + 1 } return(list(X=X, Y=Y, Z=Z, K=K)) } # test Mv1 <- mandelbrot_vectorized_1() all.equal(Mv1$K == M$K, Mv1$Z == M$Z) # TRUE
参考: https://www.r-bloggers.com/the-mandelbrot-set-in-r/
ベクトル化したRコード(2)
mandelbrot_vectorized2 <- function(xmin=-2.0, xmax=0.5, ymin=-1.25, ymax=1.25, nx=500, ny=500, max.iter=128) { X <- seq(xmin, xmax, length.out=nx) Y <- seq(ymin, ymax, length.out=ny) Z <- rep(0.0 + 0.0i, nx*ny) K <- rep(0.0, nx*ny) for (i in 1:(nx*ny)) # 2つのforループをまとめる { tmp <- mandelbrot_1(X[i%%nx + ifelse(i%%nx,0,nx)] + 1i*Y[i%/%nx + ifelse(i%%nx,1,0)], max.iter) Z[i] <- tmp$z K[i] <- tmp$k } return(list(X=X, Y=Y, Z=matrix(Z, nrow=nx, ncol=ny), K=matrix(K, nrow=nx, ncol=ny))) } # test Mv2 <- mandelbrot_vectorized2() all.equal(Mv2$K == M$K, Mv2$Z == M$Z) # TRUE
doMC (doSNOW)、foreachを使った並列化
library(foreach) library(doMC) registerDoMC(cores=as.integer(system("nproc", intern=TRUE))-1) mandelbrot_parallel <- function(xmin=-2.0, xmax=0.5, ymin=-1.25, ymax=1.25, nx=500, ny=500, max.iter=128) { X <- seq(xmin, xmax, length.out=nx) Y <- seq(ymin, ymax, length.out=ny) tmp <- foreach(i=1:(nx*ny), .combine=cbind) %dopar% { mandelbrot_1( X[i%%nx + ifelse(i%%nx,0,nx)] + 1i*Y[i%/%nx + ifelse(i%%nx,1,0)], max.iter) } Z = as.vector(tmp[1,]) K = as.vector(tmp[2,]) return(list(X=X, Y=Y, Z=matrix(Z, nrow=nx, ncol=ny), K=matrix(K, nrow=nx, ncol=ny))) } # test Mp <- mandelbrot_parallel() all.equal(Mp$K == M$K, Mp$Z == M$Z) # TRUE
ここまでの関数のパフォーマンス
ノートPCだと遅いので結局計測していない。特に並列化したやつはいつまで経っても計算が終わらないので最適化失敗である。
R | ベクトル化(2次元) | ベクトル化(1次元) | 並列化 | |
時間 | sec | sec | sec | sec |
func <- c(mandelbrot_naive, mandelbrot_vectorized1, mandelbrot_vectorized2, mandelbrot_parallel) for (f in func) { time <- foreach(i=1:10, .combine=c) %do% { system.time(f()) } print(mean(time)) }