【发布时间】:2014-02-19 19:38:35
【问题描述】:
我的问题:
我在 64 位 Windows7 PC 上使用 R.3.0.1 和 RStudio 0.97.551,并且我已经开始使用 Rcpp 将功能外包给 C/C++。该函数可以编译,但在 R 函数中对其进行评估会产生运行时错误。我不知道为什么以及如何解决这个问题。
详情
下面是我的 cpp 文件...假设它叫做“vector.cpp”
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector l5(double z, double k, double s, double K, double theta, double x, double h, NumericVector m){
int n = m.size();
NumericVector a(n);
NumericVector bu(n);
NumericVector b(n);
NumericVector c(n);
for(int i=0; i<n+1; i++){
a[i] = pow(z,m[i]) * (pow((x*pow(h,m[i])/K), theta) - 1) * (K/theta);
for (int j=0; j<i; j++){
bu[i] += pow(z,j) * (1 - z) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, j)/K), theta) - 1) * (K/theta)));
}
b[i] = k *bu[i];
c[i] = k * pow(z, m[i]) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, m[i])/K), theta) - 1) * (K/theta)));
}
return wrap(a-b-c);
}
我使用命令在 R(或 RStudio)中编译
sourceCpp(<path to file>/vector.cpp)
它编译 - 到目前为止一切顺利。但是,当我继续在其他 R 函数中使用函数 l5 时,它经常导致 R 崩溃(在 RStudio 和普通 R GUI 中)。事实上,评估它本身也不再稳定。为了重现这一点,例如尝试多次评估 l6
l6 <- function(zs, ks, ss, Ks, thetas, xs, hs){
z=zs
k=ks
s=ss
K=Ks
theta=thetas
x=xs
h=hs
m=0:30
res <- l5(z, k, s, K, theta,x, h, m)
return(res)
}
然后运行
l6(0.9, 0.1, 67, 40, 0.5, 44, 1.06)
具体来说,它会产生以下运行时错误
This application has requested Runtime to terminate it in an unusual way.
那么我的函数出了什么问题?
解决方案
正如 Dirk 下面所建议的,在 for 循环中存在一个基本错误,其中 i 从 0 运行到 n,因此有 n+1 个元素,但我只初始化了长度为 n 的向量。为了避免这个错误,我现在使用迭代器实现了这个函数
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector l5(double z, double k, double s, double K, double theta, double x, double h, NumericVector m){
int n = m.size();
NumericVector a(n);
NumericVector bu(n);
NumericVector b(n);
NumericVector c(n);
for(NumericVector::iterator i = m.begin(); i != m.end(); ++i){
a[*i] = pow(z, m[*i]) * (pow((x*pow(h, m[*i])/K), theta) - 1) * (K/theta);
for(int j=0; j<*i; j++){
bu[*i] += pow(z,j) * (1 - z) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, j)/K), theta) - 1) * (K/theta)));
}
b[*i] = k *bu[*i];
c[*i] = k * pow(z, m[*i]) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, m[*i])/K), theta) - 1) * (K/theta)));
}
return wrap(a-b-c);
}
再次感谢!
【问题讨论】:
-
迭代器解决方案是一个有趣的想法。我会改为
for (int i=0; i<n; i++),这可能是所有 C 和 C++ 中使用最广泛的行。