Michael Braun
2007-Jan-04 21:09 UTC
[Rd] Parameter changes and segfault when calling C code through .Call
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
braunm at mit.edu
(617) 253-3436
Seth Falcon
2007-Jan-04 21:31 UTC
[Rd] Parameter changes and segfault when calling C code through .Call
"Michael Braun" <braunm at MIT.EDU> writes:> 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 .CallArguments sent via the .Call interface are not copied. In almost all cases, they should be treated as read-only data. You can copy a given SEXP using duplicate. Something along these lines (untested): SEXP sigma_copy; PROTECT(sigma_copy = duplicate(sigmax)); But it seems pretty clear that sigmaView is holding a pointer to sigmaAr which in turn points to the data in sigmax. + seth