R cpp против C - выброс - PullRequest
Новая
       12

R cpp против C - выброс

2 голосов
/ 24 марта 2020

У меня есть общий вопрос с конкретным примером c. Я написал функцию для вычисления дисперсии матрицы по столбцам в C (используя интерфейс .Call) и C ++ (используя интерфейс R cpp). Глядя на следующие тесты, я задаюсь вопросом: 

> microbenchmark(times = 1000,
+                colVar(AB), # .Call Interface
+                colV(AB, ncol(AB), nrow(AB)), #Rcpp
+                apply(AB, 2, var)) #R
Unit: milliseconds
                         expr       min        lq      mean    median        uq        max neval
                   colVar(AB)  3.245000  3.350793  3.474891  3.433126  3.543796   5.110652  1000
 colV(AB, ncol(AB), nrow(AB))  4.064942  4.408336 10.215952  5.934169  6.383477  99.651530  1000
            apply(AB, 2, var) 28.260730 30.740058 46.674155 31.464449 33.586160 129.343892  1000
>

В распределении и значении функции C и C ++ работают примерно одинаково, однако, когда дело доходит до максимального значения, существует огромная разница. Кто-нибудь может объяснить мне, почему? Это особенно интересно, поскольку я пытаюсь изучать C / C ++, а также потому, что я хочу писать более сложные функции на C / C ++, где это может иметь значение. AB - матрица с размером 1000 x 1000, созданная со значениями 1 000 000 rnorm (). Ниже вы найдете коды для моих функций C и R cpp:

C (R-уровень): 

colVar <- function(x){
  .Call("colV", x, ncol(x), nrow(x))
}

C (C -Level ): 

#include <R.h>
#include <Rinternals.h>
#include <math.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 xSq[(d)];
    SEXP result;
    PROTECT(result = allocVector(REALSXP, (*nc)));
    memset(REAL(result), 0, (*nc) * sizeof(double));
    double *colVar = REAL(result);
    int fr = ((*nr) - 1);


    for(z = 0; z < (d); z++){
        xSq[z] = pow(x[z], 2);
    }

    for(i = 0; i < (*nc); i++){
        double colMean = 0;
        double xSm = 0;
        double colMsq = 0;
        for(j = 0; j < (*nr); j++){
            colMean += ((x[(j + ((*nr) * i)) ]) / (*nr));
            xSm += (xSq[(j + (*nr * i))]);
        }
        colMsq = (*nr) * (pow(colMean, 2));
        colVar[i] = ((xSm - colMsq) / fr);
    }
    UNPROTECT(1);
    return(result);
}

И функция R cpp: 

cppFunction(plugins = "unwindProtect",'NumericVector colV(NumericVector y, int n, int r){
            int nc = n;
            NumericVector x = y;
            int nr = r;
            int d = n * r;
            int i, j, z;

            // NumericVector colMean (nc);
            NumericVector xSq (d);
            // NumericVector colMsq (nc);
            // NumericVector xSm (nc);

            NumericVector colVar (nc);

            int fr = ((nr) - 1);


            for(z = 0; z < (d); z++){
               xSq[z] = x[z] * x[z];
            }

            for(i = 0; i < (nc); i++){
                double colMean = 0;
                double xSm = 0;
                double colMsq = 0;
                for(j = 0; j < (nr); j++){
                    colMean += ((x[(j + ((nr) * i)) ]) / (nr));
                    xSm += (xSq[(j + (nr * i))]);
                }
                colMsq = (nr) * (colMean * colMean);
                colVar[i] = ((xSm - colMsq) / fr);
            }
            return colVar;
            }')

Я закомментировал материал в функции C ++, чтобы сделать его максимально похожим на C функция. Если кто-нибудь из вас может помочь мне с моим вопросом, я был бы очень благодарен.

...