2013-08-28 16 views
5

Ho scritto un pezzo di codice C che dichiara una matrice quadrata di dimensione 4x4. Quindi campiona da una funzione di campionamento denominata rgig nel pacchetto GeneralizedHyperbolic in R. Inceppa la matrice utilizzando una libreria gsl da gnu e sputa il risultato. Questo è un esercizio di chiamare R da C.in esecuzione R incorporato in C

#include <stdio.h> 
#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
#include <stddef.h> 

// for gsl 
#include <gsl/gsl_machine.h> 
#include <gsl/gsl_rng.h> 
#include <gsl/gsl_randist.h> 
#include <gsl/gsl_cdf.h> 
#include <gsl/gsl_cblas.h> 
#include <gsl/gsl_sf_gamma.h> 
#include <gsl/gsl_vector.h> 
#include <gsl/gsl_matrix.h> 
#include <gsl/gsl_blas.h> 
#include <gsl/gsl_linalg.h> 

// for R embedding in C 
#include <Rinternals.h> 
#include <Rdefines.h> 
#include <Rembedded.h> 
#include <R_ext/Parse.h> 

void gsl_square_matrix_inverse (gsl_matrix *, gsl_matrix *, int); 
SEXP get_rInvGauss(void); 

int main(void) 
{ 
    // Define the dimension n of the matrix 
    // and the signum s (for LU decomposition) 
    int s, i, j, n = 4; 
    // Define all the used matrices 
    gsl_matrix * m = gsl_matrix_alloc (n, n); 
    gsl_matrix * inverse = gsl_matrix_alloc (n, n); 

    // R embedding in C 
    char *localArgs[] = {"R", "--no-save","--silent"}; 
    SEXP rInvGauss; 

    // init R embedding 
    Rf_initEmbeddedR(3, localArgs); 

    printf("\n Printing matrix m before set. size %d by %d... \n", n, n); 
    for(i=0; i<n; i++){ 
     printf("\n"); 
     for(j=0; j<n; j++){ 
      printf(" %f ", gsl_matrix_get(m, i, j)); 
     } 
    } 

    // set diagonal elements of matrix m from Inverse Gaussian Random samples 
    for(i=0; i<n; i++){ 
     rInvGauss = get_rInvGauss(); 
     gsl_matrix_set(m, i, i, *REAL(rInvGauss)); 
    } 

    Rf_endEmbeddedR(0); // end the R embedding in C 

    printf("\n Printing matrix m ..... \n"); 
    for(i=0; i<n; i++){ 
     printf("\n"); 
     for(j=0; j<n; j++){ 
      printf(" %f ", gsl_matrix_get(m, i, j)); 
     } 
    } 

    // inverse of matrix m 
    gsl_square_matrix_inverse (m, inverse, n); 

    printf("\n Printing inverse of matrix m ..... \n"); 
    for(i=0; i<n; i++){ 
     printf("\n"); 
     for(j=0; j<n; j++){ 
      printf(" %f", gsl_matrix_get(inverse, i, j)); 
     } 
    } 

    return 0; 
} 


SEXP get_rInvGauss(void) { 
    SEXP e, s, t, tmp, result; 
    int errorOccurred, n=1; 
    double chi=5, psi=4, lambda=0.5; 

    // create and evaluate 'require(GeneralizedHyperbolic)' 
    PROTECT(e = lang2(install("require"), mkString("GeneralizedHyperbolic"))); 
    R_tryEval(e, R_GlobalEnv, &errorOccurred); 
    if (errorOccurred) { 
     // handle error 
     printf("\n Error loading library GeneralizedHyperbolic:"); 
    } 
    UNPROTECT(1); 


    // Create the R expressions using a paired list 
    // rgig(n = 1, chi = 5, psi = 4, lambda = 0.5) with the R API. 
    PROTECT(t = s = allocVector(LANGSXP, 5)); 
    // could also be done by: PROTECT(t = s = allocList(5)); SET_TYPEOF(s, LANGSXP); 

    tmp = findFun(install("rgig"), R_GlobalEnv); 
    if(tmp == R_NilValue) { 
     printf("No definition for function rgig.\n"); 
     UNPROTECT(1); 
     exit(1); 
     } 
    SETCAR(t, tmp); t = CDR(t); 
    SETCAR(t, ScalarInteger(n)); SET_TAG(t, install("n")); t= CDR(t); 
    SETCAR(t, ScalarReal(chi)); SET_TAG(t, install("chi")); t= CDR(t); 
    SETCAR(t, ScalarReal(psi)); SET_TAG(t, install("psi")); t= CDR(t); 
    SETCAR(t, ScalarReal(lambda)); SET_TAG(t, install("lambda")); t= CDR(t); 
    PROTECT(result = R_tryEval(VECTOR_ELT(s, 0), R_GlobalEnv, NULL)); 
    UNPROTECT(2); 

    return(result); 

} 

void gsl_square_matrix_inverse (gsl_matrix *m, gsl_matrix *inverse, int n){ 

    int s, i, j; 
    gsl_permutation * perm = gsl_permutation_alloc (n); 
    // Make LU decomposition of matrix m 
    gsl_linalg_LU_decomp (m, perm, &s); 
    // Invert the matrix m 
    gsl_linalg_LU_invert (m, perm, inverse); 

} 

ho compilato il codice utilizzando:

R CMD SHLIB -lgsl -lgslcblas embedR_matinv.c 

con uscita:

gcc -arch x86_64 -std=gnu99 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG -I/usr/local/include -fPIC -g -O2 -c embedR_matinv.c -o embedR_matinv.o 
gcc -arch x86_64 -std=gnu99 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/lib -o embedR_matinv.so embedR_matinv.o -lgsl -lgslcblas -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation 

Quando presento:

R CMD embedR_matinv 

dà l'errore:

/Library/Frameworks/R.framework/Resources/bin/Rcmd: line 62: exec: embedR_matinv: not found 

Cosa sto sbagliando?

Ho anche cambiato il main()-test() e ha fatto un oggetto condiviso come

R CMD SHLIB -lgsl -lgslcblas embedR_matinv.c -o embedR_matinv 

con l'uscita:

gcc -arch x86_64 -std=gnu99 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG -I/usr/local/include -fPIC -g -O2 -c embedR_matinv.c -o embedR_matinv.o 

Se faccio un dyn.load("embedR_matinv.so") in R Studio, il codice viene eseguito senza alcuna interruzione cioè si blocca!

Qualche suggerimento su cosa c'è di sbagliato nel codice?

+0

Cosa viene generato quando si compila 'embedR_matinv.c'? Sei sicuro di dover emettere 'R CMD embedR_matinv', non' R CMD embedR_matinv.o' o 'R CMD embedR_matinv.so'? Forse dovresti specificare un percorso completo per il file generato? – mangusta

+0

Dove nel codice si blocca? Odora come se fosse bloccato nell'interprete R nel codice R incorporato. –

risposta

0

Stai usando un sistema Mac OSX?

provare questo:

sudo cp /Library/Frameworks/R.framework/Resources/library/Rserve/libs/x86_64/Rserve-bin.so\ /Library/Frameworks/R.framework/Resources/bin/Rserve