<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/make-your-r-simulation-models-20-times-faster/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.seascapemodels.org/rstats/2017/02/26/speeding-up-sims.html"> Bluecology blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Make your R simulation models 20 times faster</h2>
<p>R can be frustratingly slow if you use its loops. However, you can speed<br>
it up significantly (e.g. 20 times!) using the <code>Rcpp</code> package. That<br>
could turn a day long simulation into an hour long simulation.</p>
<p>I had heard years ago about the <code>Rcpp</code> package. <code>Rcpp</code> lets you write<br>
and use C ++ functions in R. However, I had never bothered to learn how<br>
to write C ++. Instead if my simulation model got too slow, I just<br>
redefined the problem (e.g. by using smaller spatial domain) so I could<br>
continue with R.</p>
<div>
<img alt="" src="https://i2.wp.com/www.seascapemodels.org/Images/speeding-up-sims/popn-tseries-1.png?w=456" data-recalc-dims="1">
<p> Output of one of the population model runs showing solutions from an R function and an Rcpp function. The C++ version gives an identical results and was up to 20 times faster. </p>
</div>
<p>I persisted with R, rather than use another language, because of its<br>
powerful graphics and the convenience of using a functional language<br>
like R to perform sensitivity analyses. More on this later.</p>
<p>The other day I was browsing Wickhams <a target="_blank" rel="nofollow" href="http://adv-r.had.co.nz/">Advanced<br>
R</a> book and realised it is actually pretty easy<br>
to write basic C++ loops.</p>
<p>Then I wondered if it would still be faster if you had to make repeated<br>
calls to the same C++ function, for instance if you wanted to run a<br>
sensitivity analysis, varying some model parameters. I like to use R for<br>
this task because the <code>purrr</code> package makes it incredibly easy to run<br>
arbitrary combinations of parameters through a function. Then it is<br>
straightforward to summarize and plot the results with <code>ggplot2</code>.</p>
<p>Turns out you can get a massive improvement, even for repeated calls to<br>
the same function. Here is a test.</p>
<h3>A population model in R</h3>
<p>First up, let’s write a simple model for simulating population size over<br>
time, according to the logistic function. The below function just takes<br>
your standard <code>r</code> (intrinsic population growth) and <code>k</code> (carrying<br>
capacity) parameters and simulates population size starting at <code>yinit</code><br>
over <code>t</code> years.</p>
<p>Further, to I have included a stochastic process, whose variation is<br>
controlled by <code>thetasd</code>, to illustrate <code>Rcpp</code> random number generator.</p>
<div>
<pre><code>logmodr &lt;- function(t, yinit, r, k, thetasd){
    y &lt;- numeric(t)
    y[1] &lt;- yinit
    theta &lt;- rnorm(t, 0, thetasd)
    for(i in 2:t){
        y[i] &lt;- y[i-1]*(r - r*(y[i-1]/k)) * exp(theta[i])
    }
    return(y)
}
</code></pre>
</div>
<p>Note that I also ran these models without the stochastic component. The<br>
speedup was even greater when you compared C++ to R without the<br>
stochastic step (about 20 times).</p>
<h3>A population model in C++</h3>
<p>Now let’s write the equivalent C++ function. You will need to install<br>
the <code>Rcpp</code> package. Note that it has some other software dependencies,<br>
so I recommend you read the guide on<br>
<a target="_blank" rel="nofollow" href="https://cran.r-project.org/web/packages/Rcpp/index.html">CRAN</a>.</p>
<p>We write the function definition as a string and pass it to<br>
<code>cppFunction</code> from <code>Rcpp</code>:</p>
<div>
<pre><code>library(Rcpp)
    cppFunction(&quot;NumericVector logmodc(int t, double yinit, double r,
double k, double thetasd){
            NumericVector y(t);
            y[0] = yinit;
      NumericVector theta = rnorm(t, 0, thetasd);
            for (int i = 1; i &lt; t; ++i){
                y[i] = y[i-1]*(r - r*(y[i-1] / k)) * exp(theta[i]);
                }
            return y;
    }
    &quot;)
</code></pre>
</div>
<p>Hopefully you can understand this, even if you are not familiar with<br>
C++. The syntax is reasonably similar to R. If you learned to program in<br>
R you may notice a few discrepencies.</p>
<p>First, C++ requires that you specify the type of each variable when its<br>
created. You can’t just create new variables without assigning them a<br>
type, and you can’t just change the type. This makes C++ more efficient<br>
than R, because the computer knows exactly how much memory to allocate a<br>
variable and doesn’t have to watch for changes.</p>
<p>Second, notice I start the iterator at time-step <code>1</code>, whereas in the R<br>
code we started at time-step <code>2</code>. In C++ vectors are indexed starting at<br>
<code>0</code>.</p>
<p>Finally, don’t forget to end lines with <code>;</code> (you can use <code>;</code> to end<br>
lines in R, but it is not essential).</p>
<h3>Running a single simulation</h3>
<p>First up, we need to define the model parameters:</p>
<div>
<pre><code>t &lt;- 100
yinit &lt;- 1
k &lt;- 20
thetasd &lt;- 0.1
r &lt;- 0.2
</code></pre>
</div>
<p>Now we can run our model. I am just going to plug the models straight<br>
into <code>microbenchmark</code>, so I can compare their times.</p>
<div>
<pre><code>library(microbenchmark)
mb1 &lt;- microbenchmark(
    logmodc(t, yinit, 1.4, k, thetasd),
    logmodr(t, yinit, 1.4, k, thetasd)
)
mb1

## Unit: microseconds
##                                expr     min       lq      mean   median
##  logmodc(t, yinit, 1.4, k, thetasd)  10.051  11.1100  12.70373  11.7415
##  logmodr(t, yinit, 1.4, k, thetasd) 179.053 198.8315 251.48889 224.3450
##        uq      max neval cld
##   12.8825   67.801   100  a
##  296.1400 1098.436   100   b
</code></pre>
</div>
<p>So the C++ version is 19 times faster.</p>
<h3>Running multiple simulations</h3>
<p>So C++ is faster for a single call to a function (that contains a loop).<br>
No surprises there. What if we want to make repeated calls to the same<br>
function, is C++ still faster than R? We might want to make repeated<br>
calls if we want to run different values of <code>r</code> through our model to do<br>
a sensitivty analysis.</p>
<p>We could increase the scope of the C++ code to include a loop over<br>
different values of <code>r</code>. However, then we would lose some of the<br>
convenience of R, which is good at manipulating data. We also wouldn’t<br>
be able to use <code>purrr</code> package to make sensitivity analysis easy.</p>
<p>First, up let’s create a sequence of <code>r</code> values:</p>
<div>
<pre><code>rseq &lt;- seq(1.1, 2.2, length.out = 10)
</code></pre>
</div>
<p>Now we can run our two models. I will use <code>purrr::map</code> (the <code>::</code> just<br>
means <code>map</code> is in the package <code>purrr</code> and avoids another call to<br>
<code>library()</code>). We will also use <code>set.seed()</code> to make sure both algorithms<br>
generate the same series of random numbers, that way we can check<br>
whether the results are identical.</p>
<div>
<pre><code>set.seed(42)
yc &lt;- purrr::map(rseq, ~logmodc(t, yinit, .x, k, thetasd))
set.seed(42)
yr &lt;- purrr::map(rseq, ~logmodr(t, yinit, .x, k, thetasd))
</code></pre>
</div>
<p><code>map</code> iteratively steps through <code>rseq</code> replacing the <code>.x</code> in the<br>
function call with each value of <code>r</code> in turn. Note that we also have to<br>
turn the function call into a formula (with <code>~</code>) to iterate in this way.</p>
<p><code>map</code> returns a list, where each element is a time-series of population<br>
sizes for a given value of <code>r</code>.</p>
<p>Let’s plot the result, for the second value of <code>r</code>:</p>
<div>
<pre><code>plot(yr[[2]], type = &quot;l&quot;, col = &quot;DarkBlue&quot;, lwd = 2)
points(yc[[2]], pch = 16, col = &quot;Tomato&quot;, cex = 0.8)
legend('topleft', legend = c(&quot;R solution&quot;,&quot;C solution&quot;),
       pch = 16, col = c(&quot;DarkBlue&quot;, &quot;Tomato&quot;))
</code></pre>
</div>
<p><img alt="" src="https://i2.wp.com/www.seascapemodels.org/Images/speeding-up-sims/popn-tseries-1.png?w=456" data-recalc-dims="1"></p>
<p>They look identical, excellent.</p>
<p>Now, let’s compare the time. Remember I had wondered whether repeated<br>
calls to a C++ function might lose some of the performance gain:</p>
<div>
<pre><code>mb2 &lt;- microbenchmark(
    purrr::map(rseq, ~logmodc(t, yinit, .x, k, thetasd)),
    purrr::map(rseq, ~logmodr(t, yinit, .x, k, thetasd))
)
mb2

## Unit: microseconds
##                                                  expr      min        lq
##  purrr::map(rseq, ~logmodc(t, yinit, .x, k, thetasd))  151.421  166.4165
##  purrr::map(rseq, ~logmodr(t, yinit, .x, k, thetasd)) 1890.549 2047.6310
##       mean    median       uq      max neval cld
##   199.9101  179.5795  221.885  371.192   100  a
##  2543.3459 2233.7455 2534.350 9173.440   100   b
</code></pre>
</div>
<p>Turns out we still gain a 12 times improvement when using C++.</p>
<p>I don’t believe I have been wasting so many hours waiting for<br>
simulations to run all these years. Learning a bit of C++ is well worth<br>
the investment.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.seascapemodels.org/rstats/2017/02/26/speeding-up-sims.html"> Bluecology blog</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/spJFV3--bRs"></a></div>