今度はR言語でマンデルブロ集合を表示してみます。
マンデルブロ集合に関する記事一覧は#マンデルブロ集合をご覧ください。
1 for文を使ったRコード
C言語の書き方に近く処理フローが理解しやすいコード。ただしあまり速くない。
## マンデルブロ集合の中心座標と半径
c0 <- -0.743+0.1145i
radius <- 0.003
## 計算の反復回数上限
max.iteration <- 1600
## 画像サイズ
image.size <- 640
## カラーマップ
cols <- colorRampPalette(c("blue", "yellow", "red", "black"))(11)
### Main Routine ###
start <- Sys.time()
x <- seq(Re(c0)-radius, Re(c0)+radius, length.out=image.size)
y <- seq(Im(c0)-radius, Im(c0)+radius, length.out=image.size)
k <- matrix(0, nrow=image.size, ncol=image.size)
for (i in 1:image.size) {
for (j in 1:image.size) {
c <- complex(real=x[i], imaginary=y[j])
z <- c
for (m in 1:max.iteration) {
if (Mod(z) >= 2) break
z <- z * z + c
}
k[i, j] <- m
}
}
end <- Sys.time()
print(difftime(end, start))## Time difference of 2.139287 minsimage(x, y, k, col=cols, xlab="Re(c)", ylab="Im(c)", asp=1)
2 最適化されたRコード
for文の二重ループ(i, j)をベクトル化したRコード。
書き方変えるだけでこんなに計算時間が変わるものかと驚きます。以下の記事を参考しました。
R-bloggers: The Mandelbrot Set in R
## マンデルブロ集合の中心座標と半径
c0 <- -0.743+0.1145i
radius <- 0.003
## 計算の反復回数上限
max.iteration <- 1600
## 画像サイズ
image.size <- 640
## カラーマップ
cols <- colorRampPalette(c("blue", "yellow", "red", "black"))(11)
### Main Routine ###
start <- Sys.time()
x <- seq(Re(c0)-radius, Re(c0)+radius, length.out=image.size)
y <- seq(Im(c0)-radius, Im(c0)+radius, length.out=image.size)
c <- outer(x, y*1i, '+')
z <- matrix(0.0, nrow=image.size, ncol=image.size)
k <- matrix(0.0, nrow=image.size, ncol=image.size)
for (i in 1:max.iteration) {
index <- which(Mod(z) < 2)
z[index] <- z[index] * z[index] + c[index]
k[index] <- k[index] + 1
}
end <- Sys.time()
print(difftime(end, start))## Time difference of 31.56294 secsimage(x, y, k, col=cols, asp=1, xlab="Re(c)", ylab="Im(c)")
3 furrrによる並列処理
furrrパッケージのfuture_map()を使用して並列処理を行うRコード。
本稿執筆に用いたマシン(2C4T)の場合、演算時間はfor文のRコードより遅くなってしまいました。別の高速マシン(4C8T)で実行するとfor文のRコードより早く、ベクトル化Rコードより遅いという結果でした。どうすれば演算時間が短くなるのかは別途試してみたい。
library(furrr)
## マンデルブロ集合の中心座標と半径
c0 <- -0.743+0.1145i
radius <- 0.003
## 計算の反復回数上限
max.iteration <- 1600
## 画像サイズ
image.size <- 640
## 関数定義
mandelbrot <- function(c) {
z <- c
for (n in 1:max.iteration) {
if (Mod(z) >= 2) break
z <- z * z + c
}
return(n)
}
## カラーマップ
cols <- colorRampPalette(c("blue", "yellow", "red", "black"))(11)
### Main Routine ###
start <- Sys.time()
plan(multicore)
x <- seq(Re(c0)-radius, Re(c0)+radius, length.out=image.size)
y <- seq(Im(c0)-radius, Im(c0)+radius, length.out=image.size)
c <- as.complex(outer(x, y*1i, '+'))
k <- future_map(c, mandelbrot)
k <- matrix(as.numeric(k), nrow=image.size, ncol=image.size)
end <- Sys.time()
print(difftime(end, start))## Time difference of 3.191491 minsimage(x, y, k, col=cols, asp=1, xlab="Re(c)", ylab="Im(c)")
4 イメージファイル直接操作
直接イメージファイルを操作してピクセル単位に着色を行う。
イメージファイル操作については【R言語】ピクセル単位のビットマップ操作参照のこと。
library(tidyverse)
library(magick)
## マンデルブロ集合の中心座標と半径
c0 <- -0.743+0.1145i
radius <- 0.003
## 計算の反復回数上限
max.iteration <- 1600
## 画像サイズ
image.size <- 640
## ビットマップ作成関数定義
make_bmp <- function(x, col, iter) {
x1 <- (x - min(x)) / (max(x) - min(x)) # 正規化([0,1])
rgba <- col(x1) # x1の値に応じて色に変換
rgba[x==iter, 1:3] <- 0 # 色を黒にする
rgba[x==iter, 4] <- 255
bmp <- array(as.raw(t(rgba)), dim=c(4, nrow(x), ncol(x)))
return(bmp)
}
### Main Routine ###
start <- Sys.time()
x <- seq(Re(c0)-radius, Re(c0)+radius, length.out=image.size)
y <- seq(Im(c0)-radius, Im(c0)+radius, length.out=image.size)
c <- outer(x, y*1i, '+')
z <- matrix(0.0, nrow=image.size, ncol=image.size)
k <- matrix(0.0, nrow=image.size, ncol=image.size)
for (i in 1:max.iteration) {
index <- which(Mod(z) < 2)
z[index] <- z[index] * z[index] + c[index]
k[index] <- k[index] + 1
}
end <- Sys.time()
print(difftime(end, start))## Time difference of 31.35135 secs## カラーマップ(alphaチャンネル必要)
cols <- colorRamp(c("blue", "yellow", "red", "black"), alpha=T)
make_bmp(k, cols, max.iteration) %>% # ビットマップデータ作成
image_read() %>% # magick_imageに変換
image_flip() # 縦軸反転
