关于Rcpp的问题
Questions about Rcpp
我在R中写了下面的计算。
但是，这在许多 "for loops" .
中使用很慢
我尝试使用 Rcpp 编写类似的计算。
但这是一个错误。
请更正我的代码。
# R
data <- matrix(1: 100, ncol = 5, nrow = 20)
Y <- 10
X <- Y - 1
Z <- matrix(ncol = 1, nrow = nrow(data) - X)
for (i in 1:(nrow(data) - X)){
Z[i, ] <- sum(data[i: (i + X), ])
}
> data
[,1] [,2] [,3] [,4] [,5]
[1,] 1 21 41 61 81
[2,] 2 22 42 62 82
[3,] 3 23 43 63 83
[4,] 4 24 44 64 84
[5,] 5 25 45 65 85
[6,] 6 26 46 66 86
[7,] 7 27 47 67 87
[8,] 8 28 48 68 88
[9,] 9 29 49 69 89
[10,] 10 30 50 70 90
[11,] 11 31 51 71 91
[12,] 12 32 52 72 92
[13,] 13 33 53 73 93
[14,] 14 34 54 74 94
[15,] 15 35 55 75 95
[16,] 16 36 56 76 96
[17,] 17 37 57 77 97
[18,] 18 38 58 78 98
[19,] 19 39 59 79 99
[20,] 20 40 60 80 100
> Z
[,1]
[1,] 2275
[2,] 2325
[3,] 2375
[4,] 2425
[5,] 2475
[6,] 2525
[7,] 2575
[8,] 2625
[9,] 2675
[10,] 2725
[11,] 2775
// Rcpp
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix myRcpp(NumericMatrix data, NumericVector Y) {
int X = Y - 1;
int i;
int nrow = data.nrow();
Rcpp::NumericMatrix Z();
for (i = 0; i < nrow - X; i++) {
Z[i] = Rcpp::sum(data( Range(i, (i + X)) , _ ));
}
return (Z);
}
后记
错误信息
sourceCpp("Rcpp/test.cpp")
g++ -m64 -I"C:/PROGRA~1/R/R-31~1.2/include" -DNDEBUG -I"C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include" -I"d:/RCompile/CRANpkg/extralibs64/local/include" -O2 -Wall -mtune=core2 -c test.cpp -o test.o
test.cpp: In function 'Rcpp::NumericMatrix myRcpp(Rcpp::NumericMatrix, Rcpp::NumericVector)':
test.cpp:6:17: error: cannot convert 'Rcpp::sugar::Minus_Vector_Primitive<14, true, Rcpp::Vector<14, Rcpp::PreserveStorage> >' to 'int' in initialization
test.cpp:12:8: warning: pointer to a function used in arithmetic [-Wpointer-arith]
test.cpp:12:51: error: no matching function for call to 'sum(Rcpp::Matrix<14>::Sub)'
test.cpp:12:51: note: candidates are:
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/sugar/functions/sum.h:98:32: note: template<bool NA, class T> Rcpp::sugar::Sum<13, NA, T> Rcpp::sum(const Rcpp::VectorBase<13, NA, VEC>&)
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/sugar/functions/sum.h:103:33: note: template<bool NA, class T> Rcpp::sugar::Sum<14, NA, T> Rcpp::sum(const Rcpp::VectorBase<14, NA, VEC>&)
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/sugar/functions/sum.h:108:32: note: template<bool NA, class T> Rcpp::sugar::Sum<10, NA, T> Rcpp::sum(const Rcpp::VectorBase<10, NA, VEC>&)
test.cpp:14:14: error: invalid conversion from 'Rcpp::NumericMatrix (*)() {aka Rcpp::Matrix<14> (*)()}' to 'int' [-fpermissive]
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/vector/Matrix.h:67:5: error: initializing argument 1 of 'Rcpp::Matrix<RTYPE, StoragePolicy>::Matrix(const int&) [with int RTYPE = 14, StoragePolicy = Rcpp::PreserveStorage]' [-fpermissive]
make: *** [test.o] Error 1
你有几个错误。大多数情况下，我建议解开明显造成编译器/模板错误的复杂表达式。
使用新代码，我得到：
R> sourceCpp("/tmp/so_question.cpp")
R> # R
R> data <- matrix(1: 100, ncol = 5, nrow = 20)
R> Y <- 10
R> X <- Y - 1
R> Z <- matrix(ncol = 1, nrow = nrow(data) - X)
R> for (i in 1:(nrow(data) - X)) {
+ Z[i, ] <- sum(data[i: (i + X), ])
+ }
R> #data
R> #Z
R>
R> myRcpp(data, Y)
[1] 2275 2325 2375 2425 2475 2525 2575 2625 2675 2725 2775
R>
代码如下——我将 R 和 C++ 合并到一个文件中。
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector myRcpp(NumericMatrix data, int y) {
int x = y - 1;
int i;
int n = data.nrow();
Rcpp::NumericVector z(n - x);
for (i = 0; i < n - x; i++) {
Rcpp::SubMatrix<REALSXP> sm = data( Range(i, (i + x)) , _ );
Rcpp::NumericMatrix m(sm);
double s = Rcpp::sum(m);
z[i] = s;
}
return z;
}
/*** R
# R
data <- matrix(1: 100, ncol = 5, nrow = 20)
Y <- 10
X <- Y - 1
Z <- matrix(ncol = 1, nrow = nrow(data) - X)
for (i in 1:(nrow(data) - X)) {
Z[i, ] <- sum(data[i: (i + X), ])
}
#data
#Z
myRcpp(data, Y)
*/
我通常在RcppArmadillo中做这样的子设置问题...
编辑： 为了完整起见，在 RcppArmadillo 中，循环体减少到一行，因为我们可以使用子矩阵访问器——这里更简单，我们只选择行—— - 我们对其求和（两次）并（明确地）转换为标量。
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::rowvec myArma(arma::mat data, int y) {
int x = y - 1;
int n = data.n_rows;
arma::rowvec z(n - x);
for (int i = 0; i < n - x; i++) {
z[i] = arma::as_scalar(arma::sum(arma::sum(data.rows(i, i+x))));
}
return z;
}
正如 Romain 所建议的，Rcpp 中的额外副本有一个成本，一旦我们得到更大的矩阵，我们就可以衡量：
R> library(rbenchmark)
R> data <- matrix(1:1e6, ncol=50)
R> res <- benchmark(myRcpp(data, Y), myArma(data, Y))
R> res[,1:4]
test replications elapsed relative
2 myArma(data, Y) 100 1.555 1.000
1 myRcpp(data, Y) 100 6.860 4.412
R>
正如我所见，您和答案都对每一行重复相同的计算（行总和）Y 次。所以一个向量化的解决方案，
data %>%
rowSums %>%
cumsum %>%
{. - lag(., Y)}
（magrittr 包中的%>%，dplyr 中的lag 方法）应该做你想做的，而且它甚至比Rcpp 快很多。
Unit: milliseconds
expr min lq mean median uq max neval cld
myArma(data, Y) 10.363077 11.633387 17.963248 12.233787 14.772540 131.67896 100 b
vectorized(data) 3.172276 3.284239 5.492879 3.441609 4.664644 63.96084 100 a
您提到您最终想使用 sd 而不是 sum，但即使在那种情况下缓存平方和也会节省很多时间。
data %>%
{(. - mu)^2} %>%
rowSums %>%
cumsum %>%
{. - lag(.,Y)} %>%
divide_by( Y*NCOL(data) - 1 )
以 mu 为组意味着您得到了另一个解决方案。
这样的 C++ 实现当然会比 R 更快，但我认为您无法比上面的代码更简洁。 :)
再见，
斯特凡
编辑：根据数据的结构和大小，您可能 运行 遇到从 cumsum 减去滞后 cumsum 减去抵消的问题（即，如果 cumsum 比单行总和大 16 个数量级你运行超出了有效数字）。在这种情况下，您可以用 RcppRoll 中的 roll_sum 替换该步骤。
作为解决方案，它实际上更短，我只是没有使用 RcppRoll 的经验，不能说是否存在问题或陷阱
data %>%
rowSums %>%
roll_sum(Y)
我在R中写了下面的计算。 但是，这在许多 "for loops" .
中使用很慢我尝试使用 Rcpp 编写类似的计算。 但这是一个错误。
请更正我的代码。
# R
data <- matrix(1: 100, ncol = 5, nrow = 20)
Y <- 10
X <- Y - 1
Z <- matrix(ncol = 1, nrow = nrow(data) - X)
for (i in 1:(nrow(data) - X)){
Z[i, ] <- sum(data[i: (i + X), ])
}
> data
[,1] [,2] [,3] [,4] [,5]
[1,] 1 21 41 61 81
[2,] 2 22 42 62 82
[3,] 3 23 43 63 83
[4,] 4 24 44 64 84
[5,] 5 25 45 65 85
[6,] 6 26 46 66 86
[7,] 7 27 47 67 87
[8,] 8 28 48 68 88
[9,] 9 29 49 69 89
[10,] 10 30 50 70 90
[11,] 11 31 51 71 91
[12,] 12 32 52 72 92
[13,] 13 33 53 73 93
[14,] 14 34 54 74 94
[15,] 15 35 55 75 95
[16,] 16 36 56 76 96
[17,] 17 37 57 77 97
[18,] 18 38 58 78 98
[19,] 19 39 59 79 99
[20,] 20 40 60 80 100
> Z
[,1]
[1,] 2275
[2,] 2325
[3,] 2375
[4,] 2425
[5,] 2475
[6,] 2525
[7,] 2575
[8,] 2625
[9,] 2675
[10,] 2725
[11,] 2775
// Rcpp
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix myRcpp(NumericMatrix data, NumericVector Y) {
int X = Y - 1;
int i;
int nrow = data.nrow();
Rcpp::NumericMatrix Z();
for (i = 0; i < nrow - X; i++) {
Z[i] = Rcpp::sum(data( Range(i, (i + X)) , _ ));
}
return (Z);
}
后记 错误信息
sourceCpp("Rcpp/test.cpp")
g++ -m64 -I"C:/PROGRA~1/R/R-31~1.2/include" -DNDEBUG -I"C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include" -I"d:/RCompile/CRANpkg/extralibs64/local/include" -O2 -Wall -mtune=core2 -c test.cpp -o test.o
test.cpp: In function 'Rcpp::NumericMatrix myRcpp(Rcpp::NumericMatrix, Rcpp::NumericVector)':
test.cpp:6:17: error: cannot convert 'Rcpp::sugar::Minus_Vector_Primitive<14, true, Rcpp::Vector<14, Rcpp::PreserveStorage> >' to 'int' in initialization
test.cpp:12:8: warning: pointer to a function used in arithmetic [-Wpointer-arith]
test.cpp:12:51: error: no matching function for call to 'sum(Rcpp::Matrix<14>::Sub)'
test.cpp:12:51: note: candidates are:
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/sugar/functions/sum.h:98:32: note: template<bool NA, class T> Rcpp::sugar::Sum<13, NA, T> Rcpp::sum(const Rcpp::VectorBase<13, NA, VEC>&)
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/sugar/functions/sum.h:103:33: note: template<bool NA, class T> Rcpp::sugar::Sum<14, NA, T> Rcpp::sum(const Rcpp::VectorBase<14, NA, VEC>&)
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/sugar/functions/sum.h:108:32: note: template<bool NA, class T> Rcpp::sugar::Sum<10, NA, T> Rcpp::sum(const Rcpp::VectorBase<10, NA, VEC>&)
test.cpp:14:14: error: invalid conversion from 'Rcpp::NumericMatrix (*)() {aka Rcpp::Matrix<14> (*)()}' to 'int' [-fpermissive]
C:/PROGRA~1/R/R-31~1.2/library/Rcpp/include/Rcpp/vector/Matrix.h:67:5: error: initializing argument 1 of 'Rcpp::Matrix<RTYPE, StoragePolicy>::Matrix(const int&) [with int RTYPE = 14, StoragePolicy = Rcpp::PreserveStorage]' [-fpermissive]
make: *** [test.o] Error 1
你有几个错误。大多数情况下，我建议解开明显造成编译器/模板错误的复杂表达式。
使用新代码，我得到：
R> sourceCpp("/tmp/so_question.cpp")
R> # R
R> data <- matrix(1: 100, ncol = 5, nrow = 20)
R> Y <- 10
R> X <- Y - 1
R> Z <- matrix(ncol = 1, nrow = nrow(data) - X)
R> for (i in 1:(nrow(data) - X)) {
+ Z[i, ] <- sum(data[i: (i + X), ])
+ }
R> #data
R> #Z
R>
R> myRcpp(data, Y)
[1] 2275 2325 2375 2425 2475 2525 2575 2625 2675 2725 2775
R>
代码如下——我将 R 和 C++ 合并到一个文件中。
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector myRcpp(NumericMatrix data, int y) {
int x = y - 1;
int i;
int n = data.nrow();
Rcpp::NumericVector z(n - x);
for (i = 0; i < n - x; i++) {
Rcpp::SubMatrix<REALSXP> sm = data( Range(i, (i + x)) , _ );
Rcpp::NumericMatrix m(sm);
double s = Rcpp::sum(m);
z[i] = s;
}
return z;
}
/*** R
# R
data <- matrix(1: 100, ncol = 5, nrow = 20)
Y <- 10
X <- Y - 1
Z <- matrix(ncol = 1, nrow = nrow(data) - X)
for (i in 1:(nrow(data) - X)) {
Z[i, ] <- sum(data[i: (i + X), ])
}
#data
#Z
myRcpp(data, Y)
*/
我通常在RcppArmadillo中做这样的子设置问题...
编辑： 为了完整起见，在 RcppArmadillo 中，循环体减少到一行，因为我们可以使用子矩阵访问器——这里更简单，我们只选择行—— - 我们对其求和（两次）并（明确地）转换为标量。
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::rowvec myArma(arma::mat data, int y) {
int x = y - 1;
int n = data.n_rows;
arma::rowvec z(n - x);
for (int i = 0; i < n - x; i++) {
z[i] = arma::as_scalar(arma::sum(arma::sum(data.rows(i, i+x))));
}
return z;
}
正如 Romain 所建议的，Rcpp 中的额外副本有一个成本，一旦我们得到更大的矩阵，我们就可以衡量：
R> library(rbenchmark)
R> data <- matrix(1:1e6, ncol=50)
R> res <- benchmark(myRcpp(data, Y), myArma(data, Y))
R> res[,1:4]
test replications elapsed relative
2 myArma(data, Y) 100 1.555 1.000
1 myRcpp(data, Y) 100 6.860 4.412
R>
正如我所见，您和答案都对每一行重复相同的计算（行总和）Y 次。所以一个向量化的解决方案，
data %>%
rowSums %>%
cumsum %>%
{. - lag(., Y)}
（magrittr 包中的%>%，dplyr 中的lag 方法）应该做你想做的，而且它甚至比Rcpp 快很多。
Unit: milliseconds
expr min lq mean median uq max neval cld
myArma(data, Y) 10.363077 11.633387 17.963248 12.233787 14.772540 131.67896 100 b
vectorized(data) 3.172276 3.284239 5.492879 3.441609 4.664644 63.96084 100 a
您提到您最终想使用 sd 而不是 sum，但即使在那种情况下缓存平方和也会节省很多时间。
data %>%
{(. - mu)^2} %>%
rowSums %>%
cumsum %>%
{. - lag(.,Y)} %>%
divide_by( Y*NCOL(data) - 1 )
以 mu 为组意味着您得到了另一个解决方案。
这样的 C++ 实现当然会比 R 更快，但我认为您无法比上面的代码更简洁。 :)
再见， 斯特凡
编辑：根据数据的结构和大小，您可能 运行 遇到从 cumsum 减去滞后 cumsum 减去抵消的问题（即，如果 cumsum 比单行总和大 16 个数量级你运行超出了有效数字）。在这种情况下，您可以用 RcppRoll 中的 roll_sum 替换该步骤。
作为解决方案，它实际上更短，我只是没有使用 RcppRoll 的经验，不能说是否存在问题或陷阱
data %>%
rowSums %>%
roll_sum(Y)