tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

Rでマンデルブロ集合をプロットする

マンデルブロ集合をRで計算する。並列化の練習もしてみる。

基本的にはMy Christmas Gift: Mandelbrot Set Computation In PythonHow To Quickly Compute The Mandelbrot Set In Pythonをなぞっている。GPUによる並列化はまだ出来ていない。

定義


wikipediaより

漸化式

 \begin{cases}z_{n+1}=z_{n}^{2}+c\\z_{0}=0\end{cases}

で定義される複素数列$z_{n}\in \mathbb{C}$が$n\to\infty$の極限で無限大に発散しない複素数$c$全体が作る集合を、マンデルブロ集合という(以下$\mathbb{M}$と表す)。

性質

 \exists {n} \in \mathbb{N}, |z_{n}|>2 \to z_{n} \notin \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))
}