Re: [Rd] Nested SEXP functions

From: Oleg Sklyar <osklyar_at_ebi.ac.uk>
Date: Thu, 15 Nov 2007 23:53:53 +0000

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:

## 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



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 15 Nov 2007 - 23:55:21 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 16 Nov 2007 - 10:30:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.