Restituisce un vettore dinamico da C a R

Question

Sto scrivendo un certo codice in C per essere chiamato in modo dinamico da R.

Questo codice genera un percorso di una casuale di Poisson processo fino ad un T. tempo desiderato Così ad ogni chiamata alla mia funzione C, la lunghezza del vettore restituito sarà diversa a seconda dei numeri casuali generati.

Quale struttura dati R dovrò creare? a LISTSXP? un altro?

Come posso crearlo e come posso aggiungerlo? E soprattutto come posso restituirlo a R?

Grazie per l'aiuto.

Answer 1

Sta a te decidere cosa vuoi usare come struttura temporanea, perché alla fine dovrai comunque assegnare un vettore per il risultato. Quindi qualunque cosa tu stia usando non è ciò che tornerai. Ci sono diverse possibilità:

1. uso Calloc + Realloc + Free che consente di espandere la memoria temporanea, se necessario. Una volta ottenuto il set completo, si assegna il vettore risultato e lo si restituisce.
2. se è possibile sovrastimare facilmente le dimensioni, è possibile sovra-allocare il vettore dei risultati e utilizzare SETLENGTH prima di restituirlo. Ci sono problemi con questo, però, perché il risultato rimarrà sovraallocato fino a quando non verrà duplicato in seguito.
3. È possibile utilizzare ciò che si è accennato a quale è una lista concatenata di blocchi vettoriali, ad esempio allocare e proteggere un pairlist di cui si aggiungono i vettori alla coda in base alle proprie esigenze. Alla fine assegni il vettore di restituzione e copia il contenuto dei blocchi che hai assegnato. Questo è più contorto di entrambi i precedenti.

Ognuno di essi presenta inconvenienti e vantaggi, quindi dipende davvero dalla vostra applicazione scegliere quella più adatta a voi.

Modifica: aggiunto un esempio di utilizzo dell'approccio pairlist. Consiglio comunque l'approccio Realloc dal momento che è molto più facile, ma comunque:

#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

Se siete aperti a passaggio da C a C++, si ottiene strato di Rcpp aggiunto gratuitamente. Ecco un esempio dalla pagina del pacchetto RcppExample (ancora piuttosto incomple):

#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 
} 

Come si vede, non allocazione di memoria esplicita, liberando, PROTECT/UNPROTECT ecc, e si ottiene un elenco di oggetti di prima classe R indietro.

Ci sono molti altri esempi, incluso in other SO questions such as this one.

Edit: E non ha davvero dire quello che i percorsi avrebbe fatto, ma come una semplice illustrazione, ecco codice C++ utilizzando le aggiunte Rcpp cumsum() e rpois() che si comportano proprio come fanno in 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 

E come prova, di nuovo in R, se poniamo il seme, siamo in grado di generare esattamente gli stessi numeri:

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>