2016-02-02 48 views
2

使用R,我想知道迭代評估多個輸入和輸出函數的最佳方法。我通過看到的情節動機:http://paulbourke.net/fractals/clifford/使用2個輸入/輸出來加速迭代函數

關鍵方程爲:

x_{n+1} = sin(A* y_n) + C* cos(A* x_n) 
y_{n+1} = sin(B* x_n) + D* cos(B* y_n) 

我想存儲每個迭代的結果。我猜還有比通過下面的代碼描述的循環去一個更快的方法:

#Parameters 
A <- -1.4 
B <- 1.6 
C <- 1.0 
D <- 0.7 

n_iter <- 10000000 

#Initial values 
x0 <- 0 
y0 <- 0 

#function to calculate n+1 points 
cliff <- function(x,y){ 
    c(sin(A*y) + C*cos(A*x), sin(B*x) + D*cos(B*y)) 
} 

#matrix to store results 
res_mat <- matrix(0,nrow=n_iter,ncol=2) 

#recursive loop (definitely not the fastest way to do this?) 
for (i in 2:n_iter){ 
    res_mat[i,] <- cliff(res_mat[i-1,1],res_mat[i-1,2]) 
} 

我想這實際上並不一定是一個單一的功能,但2,關於對方的輸出工作。任何洞察到更合適的方式來評估這些功能將不勝感激。我敢說我會從一些通用的編程建議中受益,而這些編程建議不一定是R特定的。

回答

4

一個選項是使用Rcpp;對於像這樣的迭代函數,每個新值都是前一次迭代值的複數函數,這通常會產生相當不錯的加速。

library(Rcpp) 
cliff.rcpp = cppFunction(" 
NumericMatrix cliff(int nIter, double A, double B, double C, double D) { 
    NumericMatrix x(nIter, 2); 
    for (int i=1; i < nIter; ++i) { 
    x(i,0) = sin(A*x(i-1,1)) + C*cos(A*x(i-1,0)); 
    x(i,1) = sin(B*x(i-1,0)) + D*cos(B*x(i-1,1)); 
    } 
    return x; 
}") 
cliff.rcpp(10, 1, 2, 3, 4) 
#    [,1]  [,2] 
# [1,] 0.0000000 0.0000000 
# [2,] 3.0000000 4.0000000 
# [3,] -3.7267800 -0.8614156 
# [4,] -3.2595913 -1.5266964 
# [5,] -3.9781665 -4.2182644 
# [6,] -1.1296464 -3.1953775 
# [7,] 1.3346977 3.2046776 
# [8,] 0.6386906 4.4230487 
# [9,] 1.4501988 -2.3914781 
# [10,] -0.3208062 0.5208984 

我們可以看到,這將返回相同的結果的代碼中的問題:

cliff.orig <- function(n_iter, A, B, C, D) { 
    #function to calculate n+1 points 
    cliff <- function(x,y){ 
    c(sin(A*y) + C*cos(A*x), sin(B*x) + D*cos(B*y)) 
    } 

    #matrix to store results 
    res_mat <- matrix(0,nrow=n_iter,ncol=2) 

    #recursive loop (definitely not the fastest way to do this?) 
    for (i in 2:n_iter){ 
    res_mat[i,] <- cliff(res_mat[i-1,1],res_mat[i-1,2]) 
    } 
    res_mat 
} 
identical(cliff.rcpp(10, 1, 2, 3, 4), cliff.orig(10, 1, 2, 3, 4)) 
# [1] TRUE 

對於原來的問題的投入,RCPP方法產生了〜50倍速度提升:

system.time(cliff.rcpp(10000000, -1.4, 1.6, 1.0, 0.7)) 
# user system elapsed 
# 0.661 0.046 0.717 
system.time(cliff.orig(10000000, -1.4, 1.6, 1.0, 0.7)) 
# user system elapsed 
# 34.591 0.245 35.040