I am experiencing some odd behavior with the .Call interface, and I am hoping
someone out there can help me with it. Below is a simple example (in that
there are R packages that do exactly what I want), but this code illustrates
the problem. Thanks in advance for any help you can provide.
Suppose I want to compute the log density of a multivariate normal distribution
using C code and the gsl library. My R program is:
dyn.load("mvnorm-logpdf.so")
x<-c(0,0,0,0,0,0)
mu<-c(0,0,0,0,0,0)
sig<-diag(6)
print(sig)
w<-.Call("R_mvnorm_logpdf",as.double(x),as.double(mu),sig,
as.integer(6))
print(sig) # sig has changed after .Call
This code takes the SEXP's that were passed from R and converts them to gsl
objects, and then calls the logpdf function:
# include <R.h>
# include <Rinternals.h>
# include <Rmath.h>
# include <gsl/gsl_matrix.h>
# include <gsl/gsl_vector.h>
# include <gsl/gsl_blas.h>
# include <gsl/gsl_linalg.h>
# include <gsl/gsl_math.h>
SEXP R_mvnorm_logpdf (SEXP xx, SEXP mux, SEXP sigmax, SEXP kx) {
int k = INTEGER(kx)[0];
double * xAr = REAL(xx);
double * muAr = REAL(mux);
double * sigmaAr = REAL(sigmax);
SEXP res;
gsl_vector_view xView = gsl_vector_view_array(xAr,k);
gsl_vector_view muView = gsl_vector_view_array(muAr,k);
gsl_matrix_view sigmaView = gsl_matrix_view_array(sigmaAr,k,k);
gsl_vector * x = &xView.vector;
gsl_vector * mu = &muView.vector;
gsl_matrix * sigma = &sigmaView.matrix;
1: double logans = gsl_MB_mvnorm_logpdf(x, mu, sigma, k); //
<-call logpdf here
PROTECT(res=allocVector(REALSXP,1));
REAL(res)[0] = logans;
UNPROTECT(1);
return(res);
}
The logpdf function is here
double gsl_MB_mvnorm_logpdf(gsl_vector * beta, gsl_vector * betaMean,
gsl_matrix * sigma, int k) {
// computes density of multivariate normal vector at vector
beta, with mean betaMean and cov sigma
double logdetSigma = 0;
double res;
double * kern;
int i, err;
// pointer to Cholesky decomp of sigma
gsl_matrix * sigmaChol = gsl_matrix_alloc(k, k); // define
matrix that will store Chol decomp
gsl_matrix_memcpy(sigmaChol, sigma);
gsl_linalg_cholesky_decomp(sigmaChol);
// compute logdet of sigma by 2*sum of log of diagomal elements
of chol decomp
for (i=0; i<k; i++) {
logdetSigma = logdetSigma +
log(gsl_matrix_get(sigmaChol,i,i));
}
logdetSigma = 2*logdetSigma;
// compute (beta-mean)' sigma^(-1) (beta-mean)
gsl_vector * x = gsl_vector_alloc(k);
2: // gsl_matrix_fprintf(stdout,sigma,"%f");
gsl_vector_memcpy(x, beta);
gsl_vector_sub(x, betaMean); // beta - betaMean
gsl_vector * y = gsl_vector_alloc(k);
gsl_vector_memcpy(y,x);
gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit,
sigmaChol, y); // y = inv(chol)*x from BLAS
gsl_blas_ddot(y,y,kern); // kern = y'y
// compute log density
res = -k*M_LN_SQRT_2PI - 0.5*(logdetSigma + *kern);
// release space
gsl_matrix_free(sigmaChol);
gsl_vector_free(x);
gsl_vector_free(y);
return(res);
} // end gsl_mvnorm_pdf
The problem is that after I make the .Call in R, the value of the sig matrix
changes (the 1 in the upper left corner changes from a 1 to a 0). Since I
don't make changes to the sigma object directly, I don't know how that could
happen. The following output is the sig matrix, before and after the .Call.
> source("pdfTest.R")
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 1
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 1 0 0 0 0
[3,] 0 0 1 0 0 0
[4,] 0 0 0 1 0 0
[5,] 0 0 0 0 1 0
[6,] 0 0 0 0 0 1
Through some debugging, I do know that it occurs somewhere in the logpdf
function (line 1: in the first func). I should also note that this function
does compute the density correctly.
Here's some other, potentially relevant information. Note the print statement
in line 2 of the logpdf function. If that print were "uncommented" and placed
somewhere *earlier* in the function, the elements of the sigma matrix are as
they should be, and the program runs with no problem (except for the change in
sig as described above). However, if the print statement is at line 2: or
later, the correct matrix elements are printed, but then the .Call crashes with
the following message:
*** caught segfault ***
address 0x6, cause 'memory not mapped'
(If I change the size of the matrix to diag(k), the 0x6 becomes 0xk). So, for
some reason, k is interpreted as a memory address. I double-checked for places
where I may have confused pointers and values, but I can't see where that would
have happened. Also, after searching through the R-devel archives, I noticed
that others have had other odd 'memory not mapped' errors in totally different
scenarios.
So I am flummoxed, and don't really know where to go from here.
Best wishes,
MB
------------------------------------------
Michael Braun
Assistant Professor of Marketing
MIT Sloan School of Management
38 Memorial Drive, E56-329
Cambridge, MA 02139
[EMAIL PROTECTED]
(617) 253-3436
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
