SSブログ

R: Rcppでマンデルブロー集合計算を高速化 [統計]

承前。やはりインタプリタではおそいので、Rcppをつかって高速化してみた。劇的に速くなった。

Mandelbrot002.png

[追記] 画像とコードを更新しました(10/19 23:20)。

コード
library(Rcpp)
library(inline)

## C++ source
mandelbrot.src <- "
Rcpp::NumericVector x(xx);
Rcpp::NumericVector y(yy);
unsigned short t = Rcpp::as<unsigned short>(threshold);
int nx = x.size(), ny = y.size();
Rcpp::NumericVector u(nx);

if (nx != ny) {
  return Rcpp::wrap(-1);
} else {
  for (int i = 0; i < nx; i++) {
    std::complex<double> z(0.0, 0.0);
    std::complex<double> c(x[i], y[i]);
    unsigned short k = 0;

    while (k < t) {
      z = z * z + c;
      if (abs(z) > 2.0) break;
      k++;
    }
    u[i] = k;
  }
  return Rcpp::wrap(u);
}
"

## Call C++
mandelbrot.c <- cfunction(signature(xx = "numeric",
                                    yy = "numeric",
                                    threshold = "numeric"),
                          mandelbrot.src, 
                          includes = "#include <complex>",
                          Rcpp = TRUE)


## R function
mandelbrot <- function(min.x, max.x, min.y, max.y,
                       px = 256,
                       py = round((max.y - min.y) / (max.x - min.x) * px),
                       threshold = 255) {
  if (min.x > max.x) {
    a <- min.x
    min.x <- max.x
    max.x <- a
  }
  if (min.y > max.y) {
    a <- min.y
    min.y <- max.y
    max.y <- a
  }
  if (px < 0) px <- -px
  if (py < 0) py <- -py

  rx <- seq(min.x, max.x, length = px)
  ry <- seq(min.y, max.y, length = py)
  xx <- rep(rx, py)
  yy <- rep(ry, each = px)

  z <- mandelbrot.c(xx, yy, threshold)
  matz <- matrix(z, nrow = px)
}

## sample parameters
px <- 768
py <- 768
threshold <- 1023
min.x <- -0.7625
max.x <- -0.7605
min.y <- -0.0860
max.y <- -0.0840

matz <- mandelbrot(min.x, max.x, min.y, max.y, px, py,
                   threshold = threshold)

png(file = "Mandelbrot%03d.png", width = 900, height = 900)
image(seq(min.x, max.x, length = nrow(matz)),
      seq(min.y, max.y, length = ncol(matz)), matz,
      col = rep(rainbow(256), (threshold + 1) / 256),
      xlab = "x", ylab = "y")
dev.off()

タグ:R Rcpp
nice!(0)  コメント(0)  トラックバック(1) 
共通テーマ:パソコン・インターネット

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 1