tag:blogger.com,1999:blog-4946490806848569840.post3947087451201608659..comments2023-04-24T16:57:22.851-07:00Comments on R snippets: Simulation speed: GNU R vs JuliaUnknownnoreply@blogger.comBlogger9125tag:blogger.com,1999:blog-4946490806848569840.post-43512434140395231522013-12-12T08:31:29.913-08:002013-12-12T08:31:29.913-08:00Please note that on modern multicore-CPU machines ...Please note that on modern multicore-CPU machines it is possible to boost R performance with “multicore” package. You may have noticed that R runs single-core, even on such modern machines. This in big part explains the dramatic difference in performance against Julia which has been designed for parallel processing and (I assume) makes use of all available cores on a host.<br />Package “multicore” enables R to utilize multiple cores in a straightforward manner. On my 8-core CPU machine I got over 5x speedup for the above cont model:<br />> REPS<-24<br />> system.time(lapply(rep(1000,REPS),cont.run, 10000, 1000, 0.005, 10.0, 0.01))<br />user system elapsed <br />19.054 0.000 19.127 <br />> library(multicore)<br />> system.time(mclapply(rep(1000,REPS),cont.run, 10000, 1000, 0.005, 10.0, 0.01))<br />user system elapsed <br />21.535 0.172 3.669<br /><br />As you can see the trick was to replace a call to “lapply” with a call to “mclapply” - its parallelized version provided by “multicore” package. Going further you may parallelize large chunks of existing code in a transparent, nearly effortless way by redefining basic functions such as “sapply” and “replicate” to use “mclapply” instead of “lapply”.<br />You can also gain some additional seconds by enabling implicit (background) JIT compilation in R:<br /><br />> library(compiler)<br />> enableJIT(3)<br />> system.time(lapply(rep(1000,REPS),cont.run, 10000, 1000, 0.005, 10.0, 0.01))<br />user system elapsed <br />17.417 0.009 17.512 <br />> system.time(mclapply(rep(1000,REPS),cont.run, 10000, 1000, 0.005, 10.0, 0.01))<br />user system elapsed <br />26.099 0.308 3.364Tomasz Olczakhttps://www.blogger.com/profile/15756806701801765727noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-44158125658733663782013-11-28T22:50:13.353-08:002013-11-28T22:50:13.353-08:00It is equivalent.
abs(r[i]) = abs(mul * ari) = abs...It is equivalent.<br />abs(r[i]) = abs(mul * ari) = abs(mul) * abs(ari) = ari<br />as ari is always positive ant mul is either 1 or -1.Bogumił Kamińskihttps://www.blogger.com/profile/06250268799809238730noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-44871084440500423472013-11-28T20:35:30.511-08:002013-11-28T20:35:30.511-08:00That is a pretty awesome improvement in performanc...That is a pretty awesome improvement in performance! Sorry to be a pedant (this does not affect the performance comparison), but:<br /><br />the R code:<br /><br />r[i] <- mul * sum(sig > tr) / (l * n)<br />tr[runif(n) < s] <- abs(r[i])<br /><br /><br />is not equivalent to the Julia code:<br /><br />ari = aris / (l * n)<br />r[i] = mul * ari<br />for j in 1:n<br /> if rand() < s<br /> tr[j] = ari<br /> end<br />end<br /><br /><br />as you don't multiply tr[j] by 'mul' or take the absolute value.Anonymoushttps://www.blogger.com/profile/16986267261219498729noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-34366447221702621562013-11-28T20:21:19.424-08:002013-11-28T20:21:19.424-08:00The line:
tr[runif(n) < s] <- abs(r[i])
ma...The line:<br /><br />tr[runif(n) < s] <- abs(r[i])<br /><br />makes it difficult.Anonymoushttps://www.blogger.com/profile/16986267261219498729noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-19255895497043583482013-11-26T12:22:31.838-08:002013-11-26T12:22:31.838-08:00I found this port on r-bloggers and all I can say ...I found this port on r-bloggers and all I can say is...blasphemer !Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-87302502257757959622013-11-25T10:22:47.910-08:002013-11-25T10:22:47.910-08:00I was thinking about it - but found no way.
Do you...I was thinking about it - but found no way.<br />Do you have any idea how to proceed?Bogumił Kamińskihttps://www.blogger.com/profile/06250268799809238730noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-6428066649787018422013-11-25T08:25:42.578-08:002013-11-25T08:25:42.578-08:00There may be a way to improvement performance by e...There may be a way to improvement performance by eliminating R's for loop with vectorizationAndrew Zhttps://www.blogger.com/profile/10108637160465346326noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-21279751422246900512013-11-24T10:27:36.770-08:002013-11-24T10:27:36.770-08:00I don't think we need to use a different versi...I don't think we need to use a different version R here. Fragmentation doesn't help and most of the improvement achieved through pqR will be ported to R very soon.<br /><br />Now, many R users have Julia installed in their machine and some are even core developper of the language, I think the Julia community is great with very smart guys and the day I'll want to switch I'll probably use Julia.<br /><br />Now comming back to this problem, one reason why many useRs stick with R despite this kind of performance is because of the rich library composed of more than 5000 packages and one of them is Rcpp. So learning Julia is good but learning Cpp and Rcpp is better investment for a useRs.<br /><br />Translating your Julia code to Cpp is straightforward<br /><br />Cpp side :<br /><br />#include <br /><br />using namespace Rcpp;<br /><br />// [[Rcpp::export]]<br />NumericVector cont_runcpp(int reps, int n, double d, double l, double s) {<br /> NumericVector tr(n);<br /> NumericVector r(reps);<br /> for (int i = 0; i < reps; i++) {<br /> double aris = 0.0;<br /> double sig = R::rnorm(0, d);<br /> int mul = 1;<br /> if (sig < 0) {<br /> sig = -sig;<br /> mul = -1;<br /> }<br /> for (int k=0; k < n; k++) {<br /> if (sig > tr(k)) {<br /> aris += 1;<br /> }<br /> }<br /> double ari = aris / (l * n);<br /> r(i) = mul * ari;<br /> for (int j=0; j < n; j++) {<br /> if (R::runif(0, 1) < s) {<br /> tr(j) = ari;<br /> }<br /> }<br /> }<br /> return r;<br />}<br /><br />Let's put this code into a file names "cont_runcpp.cpp"<br /><br />R side :<br /><br />require(Rcpp)<br />require(e1071)<br />require(rbenchmark)<br />sourceCpp("cont_run.cpp")<br /><br />cont_runRcpp <- function(burn_in = 1000, reps, n, d, l, s)<br /> kurtosis(cont_runcpp(reps, n, d, l, s)[burn_in:reps])<br /><br />R <- function() cont.run(1000, 10000, 1000, 0.005, 10.0, 0.01)<br />Rcpp <- function() cont_runRcpp(1000, 10000, 1000, 0.005, 10.0, 0.01)<br /><br />benchmark(R(),<br /> Rcpp(),<br /> columns = c("test", "elapsed", "relative", "replications"),<br /> replications = 100)<br /><br /> test elapsed relative replications<br />1 R() 81.950 3.944 100<br />2 Rcpp() 20.778 1.000 100<br /><br />This cpp can be improved I just took you julia code as is.<br /><br /><br />dickoahttps://www.blogger.com/profile/03289036717410396319noreply@blogger.comtag:blogger.com,1999:blog-4946490806848569840.post-56887782610794184182013-11-24T05:13:23.193-08:002013-11-24T05:13:23.193-08:00Few other things to watch:
https://news.ycombinat...Few other things to watch:<br /><br />https://news.ycombinator.com/item?id=6745135<br /><br />http://radfordneal.github.io/pqR/<br /><br />http://www.renjin.org/<br /><br />https://github.com/jtalbot/riposte<br /><br />Plus Revolution, TIBCO, and Tableau<br /><br />some of which look like giving a similar speedup without the more verbose Julia code, and the somewhat lonely Julia community.<br /><br />Gaius_noreply@blogger.com