[Rd] Nested SEXP functions
statmobile at gmail.com
statmobile at gmail.com
Fri Nov 16 03:16:38 CET 2007
On Thu, Nov 15, 2007 at 11:53:53PM +0000, Oleg Sklyar wrote:
> You assume that when you call any of those functions its arguments are
> protected (as it is the case when the function is called from within R).
> However, I do not see why they would be if you do not do it explicitly.
> Therefore, it is very likely that a result of one function, used in turn
> in a call to the other one, is garbage collected on the way and what you
> get can be totally wrong. The more R API calls you have within each of
> those functions, the more likely it is that the unprotected argument is
> garbage collected as the events are triggered within some of those
> calls.
>
> One should be very careful about protecting things in .Call interface.
> Now here is what I mean in a coding example:
Oleg,
Ah Ha!!! I was sure that I protected everything within each one of
the modular functions, so I couldn't understand why my results were
similar to a garbage collection issue. That is, on small examples I
would get the right result, but on large data it would keep failing
with wild numbers.
What you say makes complete sense, and is completely consistent with
the documentation. I appreciate the examples and how they emphasize
your point, and I'm even more grateful that you saved me the effort of
having to write sample code to ask my question.
After reading your advice, I just realized something. Instead of
doing this all within my main program, I could reduce the annoying
steps in my main code by writing the Gaussian kernel as another
modular function. It would essentially accomplish the same goal!
Thanks!
>
> ## R CODE #####
> z1 = .Call("interfaceCFun1", as.numeric(x)) # OK
> z2 = .Call("interfaceCFun12", as.numeric(x)) # Can be !OK
>
> SEXP f1(SEXP);
>
> // C code
> SEXP interfaceCFun(SEXP x) {
> / /x is protected, so it prints smth and returns new var to R
> return f1(x);
> }
> SEXP interfaceCFun2(SEXP x) {
> // for first call to f1 it prints smth of x and returns new var
> // to the second call, however in the second call the argument is no
> // more x, it is no more protected and can be released by GC, so it is
> // no clear what it is going to print there!
> return f1(f1(x));
>
> // this would be correct:
> SEXP y;
> PROTECT(y=f1(x));
> // now when we print in the following function we have memory
> protected
> y = f1(y)
> UNPROTECT(1);
> return y;
> }
>
> SEXP f1(SEXP x) {
> Rprintf("\n"); // can possibly trigger garbage collector
> Rprintf(REAL(x)[0]);
> return allocVector(REALSXP,2);
> }
>
> Best,
>
> Oleg
>
>
> On Thu, 2007-11-15 at 17:10 -0500, statmobile at gmail.com wrote:
> > Hey All,
> >
> > I was wondering if I could solicit a little advice. I have been
> > experiencing some quirkiness in my C code through .Call.
> > Unfortunately, my program is rather large, so I'm trying to create a
> > brief example. Yet before I do, I thought maybe a conceptual question
> > would be sufficient.
> >
> > Basically, I'm writing up a multidimensional Gaussian likelihood (with
> > spatial covariances and other enhancements). What I'm noticing is
> > that when I have nested SEXP functions I get inconsistent results.
> >
> > Here is what I am essentially trying to accomplish when looking at the
> > Gaussian kernel:
> >
> > l(beta) = (y-X*beta)^T V^{-1} (y-X*beta)
> >
> > Now in order to accomplish this, I wrote various linear algebra
> > subroutines to handle the R objects, we'll call this:
> >
> > SEXP XtimesY(SEXP X,SEXP Y); // X*Y
> > SEXP XminusY(SEXP X,SEXP Y); // X-Y
> > SEXP tX(SEXP X); // X^T
> > SEXP mycholinv(SEXP V); // Use cholesky decomposition for inverse
> >
> > Now, what I'm noticing is that if I call each routine individually
> > such as:
> >
> > pt1=XtimesY(X,beta); // X*beta
> > pt2=Xminus(Y,pt1); // Y-X*beta
> > pt3=tX(pt2); // (Y-X*beta)^T
> > pt4=mycholinv(V); //V^{-1}
> > pt5=XtimesY(pt2,pt4); // (Y-X*beta)^T V^{-1}
> > result=XtimesY(pt5,pt2); //(y-X*beta)^T V^{-1} (y-X*beta)
> >
> > Then the result is correct. But if instead I want to save some lines
> > of code, and I use:
> >
> > result=XtimesY(XtimesY(tX(XminusY(Y,XtimesY(X,beta))),mycholinv(V)),XminusY(Y,XtimesY(X,beta)))
> >
> > I don't always get the same result. Now my question is, should I
> > expect weird and ugly things to happen when nesting SEXP functions
> > such as this? Or is it even highly discouraged in general in C to do
> > something like this?
> >
> > If this should work, then I'll need to go through each one of the
> > functions and try to diagnose where the problem lies. Yet if it
> > shouldn't work, I'll stick with the first approach I have going now.
> >
> > Thanks in advance for your input!
> >
> > ______________________________________________
> > R-devel at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
>
More information about the R-devel
mailing list