Вернуть динамический вектор из C в R

Question

Я пишу некоторый код на C, который будет вызываться динамически из R.

Этот код генерирует путь случайного Пуассона процесса до желаемого времени T. Поэтому при каждом вызове моей функции C длина возвращаемого вектора будет отличаться в зависимости от генерируемых случайных чисел.

Какую структуру данных R мне нужно будет создать? a LISTSXP ? еще один?

Как я могу его создать и как к нему добавить? И особенно, как я могу вернуть его обратно в R?

Спасибо за помощь.

Answer 1

Это действительно зависит от вас, что вы хотите использовать в качестве временной структуры, потому что в конце концов вам все равно придется выделить вектор для результата. То, что вы будете использовать, не то, что вы вернете. Есть несколько возможностей:

1. используйте Calloc + Realloc + Free, что позволяет расширять временную память по мере необходимости. Получив полный набор, вы выделяете вектор результата и возвращаете его.
2. Если вы можете легко переоценить размер, вы можете перераспределить вектор результата и использовать SETLENGTH перед его возвратом. Однако есть проблемы с этим, потому что результат останется перераспределенным, пока не будет дублирован позже.
3. Вы можете использовать то, на что вы намекали, это цепочка векторных блоков, то есть выделять и защищать один pairlist, из которого вы добавляете векторы к хвосту по мере необходимости. В конце вы выделяете вектор возврата и копируете содержимое выделенных блоков. Это более запутанно, чем оба из вышеперечисленных.

Каждый из них имеет свои недостатки и преимущества, поэтому выбор приложения, наиболее подходящего для вас, действительно зависит от вашего приложения.

Редактировать: добавлен пример использования подхода с использованием pairlists. Я все еще рекомендовал бы подход Realloc, поскольку он намного проще, но тем не менее:

#define BLOCK_SIZE xxx /* some reasonable size for increments - could be adaptive, too */ 
SEXP block; /* last vector block */
SEXP root = PROTECT(list1(block = allocVector(REALSXP, BLOCK_SIZE)));
SEXP tail = root;
double *values = REAL(block);
int count = 0, total = 0;
do { /* your code to generate values - if you want to add one 
        first try to add it to the existing block, otherwise allocate new one */
    if (count == BLOCK_SIZE) { /* add a new block when needed */
        tail = SETCDR(tail, list1(block = allocVector(REALSXP, BLOCK_SIZE)));
        values = REAL(block);
        total += count;
        count = 0;
    }
    values[count++] = next_value;
} while (...);
total += count;
/* when done, we need to create the result vector */
{
    SEXP res = allocVector(REALSXP, total);
    double *res_values = REAL(res);
    while (root != R_NilValue) {
        int size = (CDR(root) == R_NilValue) ? count : BLOCK_SIZE;
        memcpy(res_values, REAL(CAR(root)), sizeof(double) * size);
        res_values += size;
        root = CDR(root);
    }
    UNPROTECT(1);
    return res;
 }

Answer 2

Если вы открыты для перехода с C на C ++, вы получаете дополнительный слой Rcpp бесплатно.Вот пример со страницы (все еще довольно неполного) пакета RcppExample:

#include <RcppClassic.h>
#include <cmath>

RcppExport SEXP newRcppVectorExample(SEXP vector) {
BEGIN_RCPP

    Rcpp::NumericVector orig(vector);                // keep a copy 
    Rcpp::NumericVector vec(orig.size());            // create vector same size

    // we could query size via
    //   int n = vec.size();
    // and loop over the vector, but using the STL is so much nicer
    // so we use a STL transform() algorithm on each element
    std::transform(orig.begin(), orig.end(), vec.begin(), ::sqrt);

    return Rcpp::List::create(Rcpp::Named( "result" ) = vec,
                              Rcpp::Named( "original" ) = orig) ;

END_RCPP
}

Как видите, нет явного выделения памяти, освобождения, PROTECT/UNPROTECT и т. Д., И вы получите список первого класса RВозврат

Есть еще много примеров, включенных в другие вопросы SO, такие как этот .

Редактировать: И вы действительно не сказали, чтовы бы это сделали, но в качестве простой иллюстрации приведем код на C ++ с использованием дополнений Rcpp cumsum() и rpois(), которые ведут себя так же, как и в R:

R> library(inline)
R> 
R> fun <- cxxfunction(signature(ns="integer", lambdas="numeric"),
+                    plugin="Rcpp",
+                    body='
+    int n = Rcpp::as<int>(ns);
+    double lambda = Rcpp::as<double>(lambdas);
+ 
+    Rcpp::RNGScope tmp;                  // make sure RNG behaves
+ 
+    Rcpp::NumericVector vec = cumsum( rpois( n, lambda ) );
+ 
+    return vec;
+ ')
R> set.seed(42)
R> fun(3, 0.3)
[1] 1 2 2
R> fun(4, 0.4)
[1] 1 1 1 2

И в качестве доказательства вернемсяв R, если мы устанавливаем начальное число, мы можем генерировать точно такие же числа:

R> set.seed(42)
R> cumsum(rpois(3, 0.3))
[1] 1 2 2
R> cumsum(rpois(4, 0.4))
[1] 1 1 1 2
R>