Hi Steve: Very, very nice. Thanks for the useful Rcpp script. I'm not surprised that a C++ version blows my humble little R function out of the water :) I noticed that the R function ran a lot more slowly when the sojourns were very long. It suggests that algorithms that entail conditional iteration are quite likely to be better off written in a compiled programming language that can communicate with R. It also shows off the capabilities of the Rcpp package.
Best regards, Dennis On Sun, Jul 31, 2011 at 8:32 PM, Steve Lianoglou <mailinglist.honey...@gmail.com> wrote: > Hi, > > I haven't been following this thread very closely, but I'm getting the > impression that the "inner loop" that's killing you folks here looks > quite simple (assuming it is the one provided below). > > How about trying to write the of this `f4` function below using the > rcpp/inline combo. The C/C++ you will need to write looks to be quite > trivial, let's change f4 to accept an x argument as a vector: > > I've defined f4 in the same way as Dennis did: > >> f4 <- function() >> { >> x <- sample(c(-1L,1L), 1) >> >> if (x >= 0 ) {return(1)} else { >> csum <- x >> len <- 1 >> while(csum < 0) { >> csum <- csum + sample(c(-1, 1), 1) >> len <- len + 1 >> } } >> len >> } > > Now, let's do some inline/c++ mojo: > > library(inline) > inc <- " > #include <stdio.h> > #include <stdlib.h> > #include <time.h> > " > > fxx <-cxxfunction(includes=inc, plugin="Rcpp", body=" > int len = 1; > int x = ((rand() % 2 ) == 0) ? 1 : -1; > int csum = x; > > while (csum < 0) { > x = ((rand() % 2 ) == 0) ? 1 : -1; > len++; > csum = csum + x; > } > > return wrap(len); > ") > > Assuming I've faithfully translated this into c++, the timings aren't > all that comparable. > > Doing 500 replicates with the pure R version: > > set.seed(123) > system.time(out <- replicate(500, f4())) > user system elapsed > 31.525 0.120 32.510 > > Doing 10,000 replicates using the fxx function doesn't even break a sweat: > > system.time(outxx <- replicate(10000, fxx())) > user system elapsed > 0.371 0.001 0.373 > > range(out) > [1] 1 1994308 > > range(outxx) > [1] 1 11909394 > > Hope I'm not too off of the mark, here. > -steve > -- > Steve Lianoglou > Graduate Student: Computational Systems Biology > | Memorial Sloan-Kettering Cancer Center > | Weill Medical College of Cornell University > Contact Info: http://cbio.mskcc.org/~lianos/contact > ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.