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?
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
Dove nel codice si blocca? Odora come se fosse bloccato nell'interprete R nel codice R incorporato. –