Le 16/08/10 15:33, Douglas Bates a écrit :
On Mon, Aug 16, 2010 at 8:10 AM, Romain Francois
<rom...@r-enthusiasts.com> wrote:
Nothing really conclusive here ...
I've added a preliminary version of SUGAR_BLOCK_3 (only for the case where
all inputs are vectors) so that one can write :
incl8<- '
inline double y_log_y(double y, double mu) {
return y ? y*log(y/mu) : 0.;
}
inline double yMu(double y, double mu) {
return 2. * (y_log_y(y, mu) + y_log_y(1.-y,1.-mu));
}
double yEta(double y, double eta, double w) {
return w * yMu(y, 1./(1. + exp(-eta)));
}
SUGAR_BLOCK_3(yeta, ::yEta)
'
code8<- '
NumericVector ans = yeta(NumericVector(yy), NumericVector(ee),
NumericVector(ww) );
return ans ;
'
I get :
$ Rscript ylogy.R
Le chargement a nécessité le package : methods
test replications elapsed user.self
1 fx0(yy, ee, ww) 500 0.953 0.947
2 fx1(yy, ee, ww) 500 2.576 2.521
3 fx2(yy, ee, ww) 500 2.887 2.843
4 fx3(yy, ee, ww) 500 2.898 2.854
5 fx4(yy, ee, ww) 500 1.273 1.241
6 fx5(yy, ee, ww) 500 1.502 1.450
7 fx6(yy, ee, ww) 500 1.408 1.356
8 fx7(yy, ee, ww) 500 1.303 1.273
9 fx8(yy, ee, ww) 500 1.452 1.389
So slightly better than fx5 but not as good as fx4
I guess there needs to be some investigation on why fx1, fx2 and fx3 don't
perform "well".
To me the interesting question is why fx0 does better than fx4 when my
guess is that it should not be as fast. Is it possible that the level
of optimization for the C compiler when compiling R is different from
the level of optimization when inline calls the C++ compiler?
Apparently not, with the attached file, which essentially is the code
that lives in family.c, I can actually get better performance than the
code in R:
$ Rscript ylogy9.R
Le chargement a nécessité le package : methods
test replications elapsed user.self
1 fx0(yy, ee, ww) 500 0.954 0.947
2 fx9(yy, ee, ww) 500 0.943 0.936
Probably just because I cache the symbols rather than let .Call retrieve
them each time.
So, is there something in family.c you don't do in fx4 ?
Is it this perhaps:
for (i = 0; i < n; i++) {
double etai = reta[i], tmp;
tmp = (etai < MTHRESH) ? DOUBLE_EPS :
((etai > THRESH) ? INVEPS : exp(etai));
rans[i] = x_d_opx(tmp);
}
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/bzoWrs : Rcpp svn revision 2000
|- http://bit.ly/b8VNE2 : Rcpp at LondonR, oct 5th
`- http://bit.ly/aAyra4 : highlight 0.2-2
library(Rcpp)
library(inline)
library(rbenchmark)
N <- 38798
set.seed(123454321)
ee <- rnorm(N)
yy <- rbinom(N, 1, 19354/N)
ww <- rep(1, N)
bb <- binomial()
fx0 <- function(yy,ee,ww) bb$dev.resids(yy, bb$linkinv(ee), ww)
fx9 <- local( {
inc9 <- '
static R_INLINE double x_d_opx(double x) {return x/(1 + x);}
static const double THRESH = 30.;
static const double MTHRESH = -30.;
static const double INVEPS = 1/DOUBLE_EPS;
#define _(X) X
inline double y_log_y(double y, double mu)
{
return (y) ? (y * log(y/mu)) : 0;
}
extern "C" SEXP binomial_dev_resids(SEXP y, SEXP mu, SEXP wt)
{
int i, n = LENGTH(y), lmu = LENGTH(mu), lwt = LENGTH(wt), nprot = 1;
SEXP ans;
double mui, yi, *rmu, *ry, *rwt, *rans;
if (!isReal(y)) {y = PROTECT(coerceVector(y, REALSXP)); nprot++;}
ry = REAL(y);
ans = PROTECT(duplicate(y));
rans = REAL(ans);
if (!isReal(mu)) {mu = PROTECT(coerceVector(mu, REALSXP)); nprot++;}
if (!isReal(wt)) {wt = PROTECT(coerceVector(wt, REALSXP)); nprot++;}
rmu = REAL(mu);
rwt = REAL(wt);
if (lmu != n && lmu != 1)
error(_("argument %s must be a numeric vector of length 1 or length
%d"),
"mu", n);
if (lwt != n && lwt != 1)
error(_("argument %s must be a numeric vector of length 1 or length
%d"),
"wt", n);
/* Written separately to avoid an optimization bug on Solaris cc */
if(lmu > 1) {
for (i = 0; i < n; i++) {
mui = rmu[i];
yi = ry[i];
rans[i] = 2 * rwt[lwt > 1 ? i : 0] *
(y_log_y(yi, mui) + y_log_y(1 - yi, 1 - mui));
}
} else {
mui = rmu[0];
for (i = 0; i < n; i++) {
yi = ry[i];
rans[i] = 2 * rwt[lwt > 1 ? i : 0] *
(y_log_y(yi, mui) + y_log_y(1 - yi, 1 - mui));
}
}
UNPROTECT(nprot);
return ans;
}
extern "C" SEXP logit_linkinv(SEXP eta) {
SEXP ans = PROTECT(duplicate(eta));
int i, n = LENGTH(eta);
double *rans = REAL(ans), *reta = REAL(eta);
if (!n || !isReal(eta))
error(_("Argument %s must be a nonempty numeric vector"), "eta");
for (i = 0; i < n; i++) {
double etai = reta[i], tmp;
tmp = (etai < MTHRESH) ? DOUBLE_EPS :
((etai > THRESH) ? INVEPS : exp(etai));
rans[i] = x_d_opx(tmp);
}
UNPROTECT(1);
return ans;
}
'
ff <- cfunction( signature(), '', inc = inc9 )
dll <- getDynLib( ff )
logit_linkinv__symbol <- dll$logit_linkinv
binomial_dev_resids__symbol <- dll$binomial_dev_resids
logit_linkinv <- function( eta ){
.Call( logit_linkinv__symbol, eta )
}
binomial_dev_resids <- function( y, mu, wt ){
.Call( binomial_dev_resids__symbol, y, mu, wt )
}
function(yy, ee, ww){
binomial_dev_resids(yy, logit_linkinv(ee), ww)
}
} )
benchmark(fx0(yy,ee,ww), fx9(yy,ee,ww),
columns=c("test", "replications", "elapsed", "user.self"),
replications = 500)
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
