В настоящее время я представляю себя C, чтобы расширить свои R-функции. Я написал функцию в C для некоторых вычислений, и они работают нормально. Но как только я пишу обертку в самом R, возникает какая-то ошибка. Рассмотрим мою C -функцию как «colV», а «ab c» - как некоторую произвольную матрицу.
Оператор (R) .Call("colV", abc, ncol(abc), nrow(abc))
работает совершенно нормально (каждый раз, независимо от того, как часто я его использую) ), тогда как
colV = function(x){
nc = ncol(x)
nr = nrow(x)
.Call("colV", x, nc, nr)
}
дает неверный результат при третьем использовании:
> colV(abc)
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> colV(abc)
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> colV(abc)
[1] -8.370087e+22 6.254796e-01 6.774042e-01 1.709462e+00 7.250386e-01
Если я снова объявлю обертку, первые два прогона будут работать нормально, а третий снова попытается выполнить то же самое результат появляется. Обратите внимание, что кроме внешнего вида на самом деле изменяется только первое значение!
Есть ли у кого-то идеи, что не так с оберткой? Как было сказано выше, если я использую только оператор .Call, всегда возвращается правильный результат:
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
Кроме того, у меня не было этой проблемы с использованием интерфейса. C, однако это не так. решение для меня из-за скорости и памяти. Заранее спасибо Erin
Редактировать: Вот код C:
#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
SEXP colV(SEXP y, SEXP n, SEXP r){
int *nc = INTEGER(n);
double *x = REAL(y);
int d = length(y);
int *nr = INTEGER(r);
int i, j, z;
//int d = nr * nc;
double colMean[(*nc)];
double xSq[(d)];
double colMsq[(*nc)];
double xSm[(*nc)];
SEXP result;
PROTECT(result = allocVector(REALSXP, (*nc)));
memset(REAL(result), 0, (*nc) * sizeof(double));
double *colVar = REAL(result);
//PROTECT(colVar = NEW_DOUBLE(nc));
int fr = ((*nr) - 1);
for(z = 0; z < (d); z++){
xSq[z] = x[z] * x[z];
}
for(i = 0; i < (*nc); i++){
for(j = 0; j < (*nr); j++){
colMean[i] += (x[(j + ((*nr) * i)) ]);
xSm[i] += (xSq[(j + (*nr * i))]);
}
colMean[i] = (colMean[i] / (*nr));
colMsq[i] = (*nr) * (colMean[i] * colMean[i]);
colVar[i] = ((xSm[i] - colMsq[i]) / fr);
}
UNPROTECT(1);
return(result);
}