<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/optimizing-code-vs-recognizing-patterns-with-3d-arrays/"></div></div><div><a data-url="http://www.r-bloggers.com/optimizing-code-vs-recognizing-patterns-with-3d-arrays/" data-counturl="http://www.r-bloggers.com/optimizing-code-vs-recognizing-patterns-with-3d-arrays/" data-count="horizontal" data-text="Optimizing Code vs Recognizing Patterns with 3D Arrays" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://gallery.rcpp.org//articles/optimizing-code-vs-recognizing-patterns-with-3d-arrays/"> Rcpp Gallery</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Intro</h2>
<p>As is the case with the majority of posts normally born into existence, there was an <a target="_blank" href="http://stackoverflow.com/questions/37684181/r-summing-up-neighboring-matrix-elements-how-to-speed-up">interesting problem</a> that arose on recently on StackOverflow. <a target="_blank" href="http://stackoverflow.com/users/6436094/steffen">Steffen</a>, a scientist at an unnamed weather service, faced an issue with the amount of computational time required by his triple loop in R. Specifically, Steffen needed to be able to sum over different continuous regions of elements within a 3D array structure with dimensions <code>2500 x 2500 x 50</code> to acquire the amount of neighbors. The issue is quite common within geography information sciences and/or spatial statistics. What follows next is a modified version of the response that I gave providing additional insight into various solutions.</p>
<h2>Problem Statement</h2>
<p>Consider an array of matrices that contain only 1 and 0 entries that are spread out over time with dimensions . Within each matrix, there are specific regions on the  plane of a fixed time, , that must be summed over to obtain their neighboring elements. Each region is constrained to a four by four tile given by points inside , where  and .</p>
<figure>
<pre><code data-lang="r"><span># Matrix Dimensions
</span><span>xdim</span> <span>=</span> <span>20</span>
<span>ydim</span> <span>=</span> <span>20</span>

<span># Number of Time Steps
</span><span>tdim</span> <span>=</span> <span>5</span>

<span># Generation of 3D Arrays
</span>
<span># Initial Neighborhoods
</span><span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span> <span>dim</span> <span>=</span> <span>c</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span>

<span># Result storing matrix
</span><span>res</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span> <span>dim</span> <span>=</span> <span>c</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span>

<span># Calculation loop over time
</span><span>for</span> <span>(</span><span>t</span> <span>in</span> <span>1</span><span>:</span><span>tdim</span><span>){</span>
  
  <span># Subset by specific rows
</span>  <span>for</span> <span>(</span><span>x</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>xdim</span><span>-2</span><span>)){</span>
    
    <span># Subset by specific columns
</span>    <span>for</span> <span>(</span><span>y</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>ydim</span><span>-2</span><span>)){</span>
      
      <span># Sum over each region within a time point
</span>      <span>res</span><span>[</span><span>x</span><span>,</span><span>y</span><span>,</span><span>t</span><span>]</span> <span>&lt;-</span> <span>sum</span><span>(</span><span>a</span><span>[(</span><span>x</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>x</span><span>+2</span><span>),(</span><span>y</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>y</span><span>+2</span><span>),</span><span>t</span><span>])</span>
      
    <span>}</span>
    
  <span>}</span>
  
<span>}</span></code></pre>
</figure>
<h3>Sample Result:</h3>
<p>Without a loss of generality, I’ve opted to downscale the problem to avoid a considerable amount of output while display a sample result. Thus, the sample result is of dimensions: <code>x = 6</code>, <code>y = 6</code>, <code>t = 2</code>. Therefore, the results of the neighboring clusters are:</p>
<figure>
<pre><code data-lang="r"><span>,</span> <span>,</span> <span>1</span>

     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span> <span>[,</span><span>4</span><span>]</span> <span>[,</span><span>5</span><span>]</span> <span>[,</span><span>6</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>    <span>0</span>   <span>10</span>   <span>10</span>    <span>0</span>    <span>0</span>
<span>[</span><span>4</span><span>,]</span>    <span>1</span>    <span>1</span>   <span>15</span>   <span>15</span>    <span>1</span>    <span>1</span>
<span>[</span><span>5</span><span>,]</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>
<span>[</span><span>6</span><span>,]</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>

<span>,</span> <span>,</span> <span>2</span>

     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span> <span>[,</span><span>4</span><span>]</span> <span>[,</span><span>5</span><span>]</span> <span>[,</span><span>6</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>    <span>0</span>   <span>10</span>   <span>10</span>    <span>0</span>    <span>0</span>
<span>[</span><span>4</span><span>,]</span>    <span>1</span>    <span>1</span>   <span>15</span>   <span>15</span>    <span>1</span>    <span>1</span>
<span>[</span><span>5</span><span>,]</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>    <span>0</span>
<span>[</span><span>6</span><span>,]</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span>    <span>1</span></code></pre>
</figure>
<p><strong>Note:</strong> Both time points, <code>t = 1</code> and <code>t = 2</code>, are the same!!! More on this interesting pattern later…</p>
<h2>Possible Solutions</h2>
<p>There are many ways to approach a problem like this. The first is to try to optimize the initial code base within R. Secondly, and one of the primary reasons for this article is, one can port over the computational portion of the code and perhaps add some parallelization component with the hope that such a setup would decrease the amount of time required to perform said computations. Thirdly, one can try to understand the underlying structure of the arrays by searching for patterns that can be exploited to create a cache of values which would be able to be reused instead of individually computed for each time.</p>
<p>Thus, there are really three different components within this post that will be addressed:</p>
<ol>
<li>Optimizing R code directly in hopes of obtaining a speedup in brute forcing the problem;</li>
<li>Porting over the R code into C++ and parallelizing the computation using more brute force;</li>
<li>Trying to determine different patterns in the data and exploit their weakness</li>
</ol>
<h2>Optimizing within R</h2>
<p>The initial problem statement has something that all users of <code>R</code> dread: a loop. However, it isn’t just one loop, it’s 3! As is known, one of the key downsides to <code>R</code> is looping. This problem has a lot of coverage – my favorite being the <a target="_blank" href="http://dirk.eddelbuettel.com/blog/2011/04/12/">straight, curly, or compiled</a> – and is primarily one of the key reasons why <code>Rcpp</code> is favored in loop heavy problems.</p>
<p>There are a few ways we can aim to optimize just the R code:</p>
<ol>
<li>Cache subsets within the loop.</li>
<li>Parallelize the time computation stage.</li>
</ol>
<h3>Original</h3>
<p>To show differences between functions, I’ll opt to declare the base computational loop given above as:</p>
<figure>
<pre><code data-lang="r"><span>cube_r</span> <span>=</span> <span>function</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>){</span>
  <span>for</span> <span>(</span><span>t</span> <span>in</span> <span>1</span><span>:</span><span>tdim</span><span>){</span>
    <span>for</span> <span>(</span><span>x</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>xdim</span><span>-2</span><span>)){</span>
      <span>for</span> <span>(</span><span>y</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>ydim</span><span>-2</span><span>)){</span>
        <span>res</span><span>[</span><span>x</span><span>,</span><span>y</span><span>,</span><span>t</span><span>]</span> <span>=</span> <span>sum</span><span>(</span><span>a</span><span>[(</span><span>x</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>x</span><span>+2</span><span>),(</span><span>y</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>y</span><span>+2</span><span>),</span><span>t</span><span>])</span>
      <span>}</span>
    <span>}</span>
  <span>}</span>
  <span>res</span>
<span>}</span></code></pre>
</figure>
<h3>Cached R</h3>
<p>The first order of business is to implement a cached value schema so that we avoid subsetting the different elements considerably.</p>
<figure>
<pre><code data-lang="r"><span>cube_r_cache</span> <span>=</span> <span>function</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>){</span>
  <span>for</span> <span>(</span><span>t</span> <span>in</span> <span>1</span><span>:</span><span>tdim</span><span>){</span>
    <span>temp_time</span> <span>=</span> <span>a</span><span>[,,</span><span>t</span><span>]</span>
    <span>for</span> <span>(</span><span>x</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>xdim</span><span>-2</span><span>)){</span>
      <span>temp_row</span> <span>=</span> <span>temp_time</span><span>[(</span><span>x</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>x</span><span>+2</span><span>),]</span>
      <span>for</span> <span>(</span><span>y</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>ydim</span><span>-2</span><span>)){</span>
        <span>res</span><span>[</span><span>x</span><span>,</span><span>y</span><span>,</span><span>t</span><span>]</span> <span>=</span> <span>sum</span><span>(</span><span>temp_row</span><span>[,(</span><span>y</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>y</span><span>+2</span><span>)])</span>
      <span>}</span>
    <span>}</span>
  <span>}</span>
  <span>res</span>
<span>}</span></code></pre>
</figure>
<h3>Parallelized R</h3>
<p>Next up, let’s implement a way to parallelize the computation over time using R’s built in <code>parallel</code> package.</p>
<figure>
<pre><code data-lang="r"><span>library</span><span>(</span><span>parallel</span><span>)</span>

<span>cube_r_parallel</span> <span>=</span> <span>function</span><span>(</span><span>a</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>,</span> <span>ncores</span><span>){</span>
  
  <span># Create a cluster
</span>  <span>cl</span> <span>=</span> <span>makeCluster</span><span>(</span><span>ncores</span><span>)</span>
  
  <span># Create a time-based computation loop
</span>  <span>computation_loop</span> <span>=</span> <span>function</span><span>(</span><span>X</span><span>,</span> <span>a</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>){</span>
    
    <span>temp_time</span> <span>=</span> <span>a</span><span>[,,</span><span>X</span><span>]</span>
    <span>res</span> <span>=</span> <span>temp_time</span>
    
    <span>for</span> <span>(</span><span>x</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>xdim</span><span>-2</span><span>)){</span>
      <span>temp_row</span> <span>=</span> <span>temp_time</span><span>[(</span><span>x</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>x</span><span>+2</span><span>),]</span>
      <span>for</span> <span>(</span><span>y</span> <span>in</span> <span>3</span><span>:</span><span>(</span><span>ydim</span><span>-2</span><span>)){</span>
        <span>res</span><span>[</span><span>x</span><span>,</span><span>y</span><span>]</span> <span>=</span> <span>sum</span><span>(</span><span>temp_row</span><span>[,(</span><span>y</span><span>-2</span><span>)</span><span>:</span><span>(</span><span>y</span><span>+2</span><span>)])</span>
      <span>}</span>
    <span>}</span>
    
    <span>res</span>
  <span>}</span>
  
  <span># Obtain the results
</span>  <span>r_parallel</span> <span>=</span> <span>parSapply</span><span>(</span><span>cl</span> <span>=</span> <span>cl</span><span>,</span>
                         <span>X</span> <span>=</span> <span>seq_len</span><span>(</span><span>tdim</span><span>),</span> 
                         <span>FUN</span> <span>=</span> <span>computation_loop</span><span>,</span> 
                         <span>a</span> <span>=</span> <span>a</span><span>,</span> <span>xdim</span> <span>=</span> <span>xdim</span><span>,</span> <span>ydim</span> <span>=</span> <span>ydim</span><span>,</span> <span>tdim</span> <span>=</span> <span>tdim</span><span>,</span> 
                         <span>simplify</span> <span>=</span> <span>&quot;array&quot;</span><span>)</span>
  
  <span># Kill the cluster
</span>  <span>stopCluster</span><span>(</span><span>cl</span><span>)</span>
  
  <span># Return
</span>  <span>r_parallel</span>
<span>}</span></code></pre>
</figure>
<h3>Benchmarking</h3>
<p>As this <strong>is</strong> <code>Rcpp</code>, benchmarks are king. Here I’ve opted to create a benchmark of the different functions after trying to optimize them within R.</p>
<figure>
<pre><code data-lang="r"><span>rtimings</span> <span>=</span> <span>microbenchmark</span><span>(</span><span>op_answer</span> <span>=</span> <span>cube_r</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>),</span>
                          <span>r_cached_subset</span> <span>=</span> <span>cube_r_cache</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>),</span>
                          <span>r_parallel1</span> <span>=</span> <span>cube_r_parallel</span><span>(</span><span>a</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>,</span> <span>ncores</span> <span>=</span> <span>1</span><span>),</span>
                          <span>r_parallel2</span> <span>=</span> <span>cube_r_parallel</span><span>(</span><span>a</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>,</span> <span>ncores</span> <span>=</span> <span>2</span><span>),</span>
                          <span>r_parallel3</span> <span>=</span> <span>cube_r_parallel</span><span>(</span><span>a</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>,</span> <span>ncores</span> <span>=</span> <span>3</span><span>),</span>
                          <span>r_parallel4</span> <span>=</span> <span>cube_r_parallel</span><span>(</span><span>a</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>,</span> <span>ncores</span> <span>=</span> <span>4</span><span>))</span></code></pre>
</figure>
<p><strong>Output:</strong></p>
<figure>
<pre><code data-lang="r"><span>Unit</span><span>:</span> <span>milliseconds</span>
            <span>expr</span>        <span>min</span>         <span>lq</span>       <span>mean</span>     <span>median</span>         <span>uq</span>        <span>max</span> <span>neval</span>
 <span>r_cached_subset</span>   <span>2.914021</span>   <span>3.380441</span>   <span>3.762666</span>   <span>3.498877</span>   <span>3.671357</span>   <span>6.584328</span>   <span>100</span>
       <span>op_answer</span>   <span>4.007368</span>   <span>4.350468</span>   <span>5.053147</span>   <span>4.630304</span>   <span>4.934394</span>   <span>9.320869</span>   <span>100</span>
     <span>r_parallel1</span> <span>143.999263</span> <span>149.833078</span> <span>152.525882</span> <span>152.473543</span> <span>154.586970</span> <span>165.265859</span>   <span>100</span>
     <span>r_parallel2</span> <span>284.159601</span> <span>291.980428</span> <span>296.283743</span> <span>295.066914</span> <span>300.471994</span> <span>312.267809</span>   <span>100</span>
     <span>r_parallel3</span> <span>426.694259</span> <span>437.278416</span> <span>441.646800</span> <span>440.875021</span> <span>445.099250</span> <span>475.146716</span>   <span>100</span>
     <span>r_parallel4</span> <span>563.040544</span> <span>582.282379</span> <span>588.159337</span> <span>585.561017</span> <span>593.944895</span> <span>644.125363</span>   <span>100</span></code></pre>
</figure>
<p>Not surprisingly, the cached subset function performed better than the original function. What was particularly surpising in this case was that the parallelization option within R took a considerably longer amount of time as additional cores were added. This was partly due to the fact that the <code>a</code> array had to be replicated out across the different <code>R</code> processes and then there was also a time lag between spawning the processes and winding them down. This may prove to be a fatal flaw of the construction of <code>cube_r_parallel</code> as a normal user may wish to make repeated calls within the same parallelized session. However, I digress as I feel like I’ve spent too much time describing ways to optimize the <code>R</code> code.</p>
<h2>Porting into C++</h2>
<p>To keep in the spirit of pure ‘brute force’, we can quickly port over the R cached version of the computational loop to <a target="_blank" href="http://arma.sourceforge.net/docs.html"><code>Armadillo</code>, C++ Linear Algebra Library,</a> and use <a target="_blank" href="http://openmp.org/wp/"><code>OpenMP</code></a> with <a target="_blank" href="https://cran.r-project.org/web/packages/RcppArmadillo/index.html"><code>RcppArmadillo</code></a>.</p>
<p>Why not use just <code>Rcpp</code> though? The rationale for using <code>RcppArmadillo</code> over <code>Rcpp</code> is principally because of the native support for multidimensional arrays via <a target="_blank" href="http://arma.sourceforge.net/docs.html#Cube">arma::cube</a>. <a target="_blank" href="http://stackoverflow.com/a/12570932/1345455">In addition, it’s a bit painful in the current verison of <code>Rcpp</code> to work with an array.</a> Hence, I’d rather face the cost of copying an object from <code>SEXP</code> to <code>arma</code> and back again than mess around with this.</p>
<p>Before we get started, it is very important to provide protection for systems that still lack <code>OpenMP</code> by default in R (**cough** OS X **cough**). Though, this option can be enabled using the steps described at <a target="_blank" href="http://thecoatlessprofessor.com/programming/openmp-in-r-on-os-x/">http://thecoatlessprofessor.com/programming/openmp-in-r-on-os-x/</a>. After saying that, we still need to protect users who do not use that option by guarding the header inclusion:</p>
<figure>
<pre><code data-lang="cpp"><span>// Protect against compilers without OpenMP
</span><span>#ifdef _OPENMP
#include &lt;omp.h>
#endif</span></code></pre>
</figure>
<p>This provides protection from the compiler throwing an error due to <code>OpenMP</code> not being available on the system. Note: In the event that the system does not have <code>OpenMP</code>, the process will be executed serially just like always.</p>
<p>With this being said, let’s look at the C++ port of the R function:</p>
<figure>
<pre><code data-lang="cpp"><span>#include &lt;RcppArmadillo.h>
</span>
<span>// [[Rcpp::depends(RcppArmadillo)]]
</span>
<span>// Add a flag to enable OpenMP at compile time
// [[Rcpp::plugins(openmp)]]
</span>
<span>// Protect against compilers without OpenMP
</span><span>#ifdef _OPENMP
#include &lt;omp.h>
#endif
</span>
<span>// Parallelized C++ equivalent
// [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>cube</span> <span>cube_cpp_parallel</span><span>(</span><span>arma</span><span>::</span><span>cube</span> <span>a</span><span>,</span> <span>arma</span><span>::</span><span>cube</span> <span>res</span><span>,</span> <span>int</span> <span>ncores</span> <span>=</span> <span>1</span><span>)</span> <span>{</span>
  
  <span>// Extract the different dimensions
</span>
  <span>// Normal Matrix dimensions
</span>  <span>unsigned</span> <span>int</span> <span>xdim</span> <span>=</span> <span>res</span><span>.</span><span>n_rows</span><span>;</span>
  
  <span>unsigned</span> <span>int</span> <span>ydim</span> <span>=</span> <span>res</span><span>.</span><span>n_cols</span><span>;</span>
  
  <span>// Depth of Array
</span>  <span>unsigned</span> <span>int</span> <span>tdim</span> <span>=</span> <span>res</span><span>.</span><span>n_slices</span><span>;</span>
  
  <span>// Added an omp pragma directive to parallelize the loop with ncores
</span>  <span>#pragma omp parallel for num_threads(ncores)
</span>  <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>t</span> <span>=</span> <span>0</span><span>;</span> <span>t</span> <span>&lt;</span> <span>tdim</span><span>;</span> <span>t</span><span>++</span><span>){</span> <span>// Note: Array's Index in C++ starts at 0 and goes to N - 1
</span>    
    <span>// Pop time `t` from `a` e.g. `a[,,t]`
</span>    <span>arma</span><span>::</span><span>mat</span> <span>temp_mat</span> <span>=</span> <span>a</span><span>.</span><span>slice</span><span>(</span><span>t</span><span>);</span>
    
    <span>// Begin region selection
</span>    <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>x</span> <span>=</span> <span>2</span><span>;</span> <span>x</span> <span>&lt;</span> <span>xdim</span><span>-</span><span>2</span><span>;</span> <span>x</span><span>++</span><span>){</span>
      
      <span>// Subset the rows
</span>      <span>arma</span><span>::</span><span>mat</span> <span>temp_row_sub</span> <span>=</span> <span>temp_mat</span><span>.</span><span>rows</span><span>(</span><span>x</span><span>-</span><span>2</span><span>,</span> <span>x</span><span>+</span><span>2</span><span>);</span>
      
      <span>// Iterate over the columns with unit accumulative sum
</span>      <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>y</span> <span>=</span> <span>2</span><span>;</span> <span>y</span> <span>&lt;</span>  <span>ydim</span><span>-</span><span>2</span><span>;</span> <span>y</span><span>++</span><span>){</span>
        <span>res</span><span>(</span><span>x</span><span>,</span><span>y</span><span>,</span><span>t</span><span>)</span> <span>=</span> <span>accu</span><span>(</span><span>temp_row_sub</span><span>.</span><span>cols</span><span>(</span><span>y</span><span>-</span><span>2</span><span>,</span><span>y</span><span>+</span><span>2</span><span>));</span>
      <span>}</span>
    <span>}</span>
  <span>}</span>
  
  <span>return</span> <span>res</span><span>;</span>
<span>}</span></code></pre>
</figure>
<p>A few things to note here:</p>
<ol>
<li><code>a</code>, <code>res</code>, <code>xdim</code>,<code>ydim</code> and <code>ncores</code> are <strong>shared</strong> across processes.</li>
<li><code>t</code> is unique to each process.</li>
<li>We protect users that cannot use <code>OpenMP</code>!</li>
</ol>
<p>To verify that this is an equivalent function, we opt to check the object equality:</p>
<figure>
<pre><code data-lang="r"><span>cpp_parallel</span> <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>ncores</span><span>)</span>

<span>all.equal</span><span>(</span><span>cpp_parallel</span><span>,</span> <span>op_answer</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>[</span><span>1</span><span>]</span> <span>TRUE</span></code></pre>
</figure>
<h3>Timings</h3>
<p>Just as before, let’s check the benchmarks to see how well we did:</p>
<figure>
<pre><code data-lang="r"><span>port_timings</span> <span>=</span> <span>microbenchmark</span><span>(</span><span>op_answer</span> <span>=</span> <span>cube_r</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>),</span>
                              <span>r_cached_subset</span> <span>=</span> <span>cube_r_cache</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>),</span>
                              <span>cpp_core1</span> <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>1</span><span>),</span>
                              <span>cpp_core2</span> <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>2</span><span>),</span>
                              <span>cpp_core3</span> <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>3</span><span>),</span>
                              <span>cpp_core4</span> <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span><span>res</span><span>,</span><span>4</span><span>))</span></code></pre>
</figure>
<p><strong>Output:</strong></p>
<figure>
<pre><code data-lang="r"><span>Unit</span><span>:</span> <span>microseconds</span>
            <span>expr</span>      <span>min</span>        <span>lq</span>      <span>mean</span>    <span>median</span>        <span>uq</span>      <span>max</span> <span>neval</span>
       <span>op_answer</span> <span>3618.381</span> <span>3865.0475</span> <span>4145.8970</span> <span>4008.6600</span> <span>4190.3205</span> <span>5965.215</span>   <span>100</span>
 <span>r_cached_subset</span> <span>2790.882</span> <span>2992.2785</span> <span>3313.1804</span> <span>3100.9775</span> <span>3296.5020</span> <span>5141.932</span>   <span>100</span>
       <span>cpp_core1</span>  <span>174.254</span>  <span>181.6835</span>  <span>194.7873</span>  <span>185.8370</span>  <span>204.2300</span>  <span>373.787</span>   <span>100</span>
       <span>cpp_core2</span>  <span>112.481</span>  <span>118.5510</span>  <span>134.7780</span>  <span>123.8365</span>  <span>145.9055</span>  <span>343.364</span>   <span>100</span>
       <span>cpp_core3</span>   <span>82.790</span>   <span>92.0120</span>  <span>124.8239</span>  <span>108.9060</span>  <span>144.7930</span>  <span>314.537</span>   <span>100</span>
       <span>cpp_core4</span>   <span>81.068</span>  <span>116.8390</span>  <span>134.6073</span>  <span>126.3345</span>  <span>149.9915</span>  <span>256.079</span>   <span>100</span></code></pre>
</figure>
<p>Wow! The C++ version of the parallelization really did wonders vs. the previous R implementation of parallelization. The speed ups are about 44x vs. original and 34x vs. the optimized R loop (using 3 cores). Note, with the addition of the 4th core, the parallelization option performs poorly as the system running the benchmarks only has four cores. Thus, one of the cores is trying to keep up with the parallelization while also having to work on operating system tasks and so on.</p>
<p>Now, this isn’t to say that there is no cap to parallelization benefits given infinite amounts of processing power. In fact, there are two laws that govern general speedups from parallelization: <a target="_blank" href="https://en.wikipedia.org/wiki/Amdahl%27s_law">Amdahl’s Law (Fixed Problem Size)</a> and <a target="_blank" href="https://en.wikipedia.org/wiki/Gustafson%27s_law">Gustafson’s Law (Scaled Speedup)</a>.</p>
<h2>Detecting Patterns</h2>
<p>Firstly, notice how the array is constructed in this case with: <code>array(0:1, dims)</code>. There seems to be some sort of pattern depending on the <code>xdim</code>, <code>ydim</code>, and <code>tdim</code>. If we can recognize the pattern in advance, the amount of computation required decreases. <strong>However, this may also impact the ability to <em>generalize</em> the algorithm to other cases outside the problem statement.</strong> Thus, the reason for this part of the post being at the terminal part of this article.</p>
<p>After some trial and error using different dimensions, different patterns become recognizable and reproducible. Most notably, we have three different cases:</p>
<ul>
<li>Case 1: If <code>xdim</code> is even, then only the rows of a matrix alternate with rows containing <strong>all</strong> 1s or 0s.</li>
<li>Case 2: If <code>xdim</code> is odd and <code>ydim</code> is even, then only the rows alternate of a matrix alternate with rows containing a <strong>combination of both</strong> 1 or 0.</li>
<li>Case 3: If <code>xdim</code> is odd and <code>ydim</code> is odd, then rows alternate as well as the matrices alternate with a <strong>combination of both</strong> 1 or 0.</li>
</ul>
<h2>Examples of Pattern Cases</h2>
<p>Let’s see the cases in action to observe the patterns. Note, the size of these arrays is small and would likely yield issues where the <code>for</code> loop given above would not work due to an out-of-bounds error.</p>
<p><strong>Case 1:</strong></p>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>2</span>
<span>ydim</span> <span>=</span> <span>3</span>
<span>tdim</span> <span>=</span> <span>2</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span> <span>dim</span> <span>=</span> <span>c</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<p><strong>Output</strong>:</p>
<figure>
<pre><code data-lang="r"><span>,</span> <span>,</span> <span>1</span>

     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>0</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>1</span>    <span>1</span>

<span>,</span> <span>,</span> <span>2</span>

     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>0</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>1</span>    <span>1</span></code></pre>
</figure>
<p><strong>Case 2:</strong></p>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>3</span>
<span>ydim</span> <span>=</span> <span>4</span>
<span>tdim</span> <span>=</span> <span>2</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span> <span>dim</span> <span>=</span> <span>c</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<p><strong>Output:</strong></p>
<figure>
<pre><code data-lang="r"><span>,</span> <span>,</span> <span>1</span>
    
     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span> <span>[,</span><span>4</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>    <span>1</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>0</span>    <span>1</span>    <span>0</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>    <span>1</span>
    
<span>,</span> <span>,</span> <span>2</span>
    
     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span> <span>[,</span><span>4</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>    <span>1</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>0</span>    <span>1</span>    <span>0</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>    <span>1</span></code></pre>
</figure>
<p><strong>Case 3:</strong></p>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>3</span>
<span>ydim</span> <span>=</span> <span>3</span>
<span>tdim</span> <span>=</span> <span>3</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span> <span>dim</span> <span>=</span> <span>c</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<p><strong>Output:</strong></p>
<figure>
<pre><code data-lang="r"><span>,</span> <span>,</span> <span>1</span>
    
     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>0</span>    <span>1</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>
    
<span>,</span> <span>,</span> <span>2</span>
    
     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>1</span>    <span>0</span>    <span>1</span>
<span>[</span><span>2</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>
<span>[</span><span>3</span><span>,]</span>    <span>1</span>    <span>0</span>    <span>1</span>
    
<span>,</span> <span>,</span> <span>3</span>
    
     <span>[,</span><span>1</span><span>]</span> <span>[,</span><span>2</span><span>]</span> <span>[,</span><span>3</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>    <span>0</span>    <span>1</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>    <span>1</span>    <span>0</span></code></pre>
</figure>
<h3>Pattern Hacking</h3>
<p>Based on the above discussion, we opt to make a bit of code that exploits this unique pattern. The language that we can write this code in is either R or C++. Though, due to the nature of this website, we opt to proceed in the later. Nevertheless, as an exercise, feel free to backport this code into R.</p>
<p>Having acknowledged that, we opt to start off trying to create code that can fill a matrix with either an <code>even</code> or an <code>odd</code> column vector. e.g.</p>
<h4>Odd Vector</h4>
<p>The <em>odd</em> vector is defined as having the initial value being given by <code>1</code> instead of <code>0</code>. This is used in heavily in both <strong>Case 2</strong> and <strong>Case 3</strong>.</p>
<figure>
<pre><code data-lang="r">     <span>[,</span><span>1</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>1</span>
<span>[</span><span>2</span><span>,]</span>    <span>0</span>
<span>[</span><span>3</span><span>,]</span>    <span>1</span>
<span>[</span><span>4</span><span>,]</span>    <span>0</span>
<span>[</span><span>5</span><span>,]</span>    <span>1</span></code></pre>
</figure>
<h4>Even</h4>
<p>In comparison to the <em>odd</em> vector, the <em>even</em> vector starts at 0. This is used principally in <strong>Case 1</strong> and then in <strong>Case 2</strong> and <strong>Case 3</strong> to alternate columns.</p>
<figure>
<pre><code data-lang="r">     <span>[,</span><span>1</span><span>]</span>
<span>[</span><span>1</span><span>,]</span>    <span>0</span>
<span>[</span><span>2</span><span>,]</span>    <span>1</span>
<span>[</span><span>3</span><span>,]</span>    <span>0</span>
<span>[</span><span>4</span><span>,]</span>    <span>1</span>
<span>[</span><span>5</span><span>,]</span>    <span>0</span></code></pre>
</figure>
<h4>Creating Alternating Vectors</h4>
<p>To obtain such an alternating vector that switches between two values, we opt to create a vector using the modulus of the vector position <code>i % 2</code>, where <code>0 % 2</code> or <code>2 % 2</code> is 0 (divisible by 2 / even) and <code>1 % 2</code> is 1 (not divisible by 2 / odd). alternating vector in this case switches between two different values.</p>
<figure>
<pre><code data-lang="cpp"><span>#include &lt;RcppArmadillo.h>
</span><span>// [[Rcpp::depends(RcppArmadillo)]]
</span>
<span>// ------- Make Alternating Column Vectors
</span>
<span>// Creates a vector with initial value 1
</span><span>arma</span><span>::</span><span>vec</span> <span>odd_vec</span><span>(</span><span>unsigned</span> <span>int</span> <span>xdim</span><span>){</span>
  
  <span>// make a temporary vector to create alternating 0-1 effect by row.
</span>  <span>arma</span><span>::</span><span>vec</span> <span>temp_vec</span><span>(</span><span>xdim</span><span>);</span>
  
  <span>// Alternating vector (anyone have a better solution? )
</span>  <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>i</span> <span>=</span> <span>0</span><span>;</span> <span>i</span> <span>&lt;</span> <span>xdim</span><span>;</span> <span>i</span><span>++</span><span>)</span> <span>{</span>
    <span>temp_vec</span><span>(</span><span>i</span><span>)</span> <span>=</span> <span>(</span><span>i</span> <span>%</span> <span>2</span> <span>?</span> <span>0</span> <span>:</span> <span>1</span><span>);</span> <span>// Ternary operator in C++, e.g. if(TRUE){1}else{0}
</span>  <span>}</span>
  
  <span>return</span> <span>temp_vec</span><span>;</span>
<span>}</span>

<span>// Creates a vector with initial value 0
// [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>vec</span> <span>even_vec</span><span>(</span><span>unsigned</span> <span>int</span> <span>xdim</span><span>){</span>
  
  <span>// make a temporary vector to create alternating 0-1 effect by row.
</span>  <span>arma</span><span>::</span><span>vec</span> <span>temp_vec</span><span>(</span><span>xdim</span><span>);</span>
  
  <span>// Alternating vector (anyone have a better solution? )
</span>  <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>i</span> <span>=</span> <span>0</span><span>;</span> <span>i</span> <span>&lt;</span> <span>xdim</span><span>;</span> <span>i</span><span>++</span><span>)</span> <span>{</span>
    <span>temp_vec</span><span>(</span><span>i</span><span>)</span> <span>=</span> <span>(</span><span>i</span> <span>%</span> <span>2</span> <span>?</span> <span>1</span> <span>:</span> <span>0</span><span>);</span> <span>// changed
</span>  <span>}</span>
  
  <span>return</span> <span>temp_vec</span><span>;</span>
<span>}</span></code></pre>
</figure>
<h4>Creating the three cases of matrix</h4>
<p>With our ability to now generate <strong>odd</strong> and <strong>even</strong> vectors by column, we now need to figureo out how to create the matrices described in each case. As mentioned above, there are three cases of matrix given as:</p>
<ol>
<li>The even,</li>
<li>The bad odd,</li>
<li>And the ugly odd.</li>
</ol>
<figure>
<pre><code data-lang="cpp"><span>// --- Requires the previous vector code snippets to be within one file inorder to run.
</span>
<span>// --- Handle the different matrix cases </span>
<span>// Case 1: xdim is even
// [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>mat</span> <span>make_even_matrix_case1</span><span>(</span><span>unsigned</span> <span>int</span> <span>xdim</span><span>,</span> <span>unsigned</span> <span>int</span> <span>ydim</span><span>){</span>
  
  <span>arma</span><span>::</span><span>mat</span> <span>temp_mat</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>);</span>
  
  <span>temp_mat</span><span>.</span><span>each_col</span><span>()</span> <span>=</span> <span>even_vec</span><span>(</span><span>xdim</span><span>);</span>
  
  <span>return</span> <span>temp_mat</span><span>;</span>
<span>}</span>

<span>// Case 2: xdim is odd and ydim is odd // [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>mat</span> <span>make_odd_matrix_case2</span><span>(</span><span>unsigned</span> <span>int</span> <span>xdim</span><span>,</span> <span>unsigned</span> <span>int</span> <span>ydim</span><span>){</span>
  
  <span>arma</span><span>::</span><span>mat</span> <span>temp_mat</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>);</span>
  
  <span>// Cache values
</span>  <span>arma</span><span>::</span><span>vec</span> <span>e_vec</span> <span>=</span> <span>even_vec</span><span>(</span><span>xdim</span><span>);</span>
  <span>arma</span><span>::</span><span>vec</span> <span>o_vec</span> <span>=</span> <span>odd_vec</span><span>(</span><span>xdim</span><span>);</span>
  
  <span>// Alternating column </span>  <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>i</span> <span>=</span> <span>0</span><span>;</span> <span>i</span> <span>&lt;</span> <span>ydim</span><span>;</span> <span>i</span><span>++</span><span>)</span> <span>{</span>
    <span>temp_mat</span><span>.</span><span>col</span><span>(</span><span>i</span><span>)</span> <span>=</span> <span>(</span><span>i</span> <span>%</span> <span>2</span> <span>?</span> <span>e_vec</span> <span>:</span> <span>o_vec</span><span>);</span> 
  <span>}</span>
  
  <span>return</span> <span>temp_mat</span><span>;</span>
<span>}</span>

<span>// Case 3: xdim is odd and ydim is even
// [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>mat</span> <span>make_odd_matrix_case3</span><span>(</span><span>unsigned</span> <span>int</span> <span>xdim</span><span>,</span> <span>unsigned</span> <span>int</span> <span>ydim</span><span>){</span>
  
  <span>arma</span><span>::</span><span>mat</span> <span>temp_mat</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>);</span>
  
  <span>// Cache values
</span>  <span>arma</span><span>::</span><span>vec</span> <span>e_vec</span> <span>=</span> <span>even_vec</span><span>(</span><span>xdim</span><span>);</span>
  <span>arma</span><span>::</span><span>vec</span> <span>o_vec</span> <span>=</span> <span>odd_vec</span><span>(</span><span>xdim</span><span>);</span>
  
  <span>// Alternating column </span>  <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>i</span> <span>=</span> <span>0</span><span>;</span> <span>i</span> <span>&lt;</span> <span>ydim</span><span>;</span> <span>i</span><span>++</span><span>)</span> <span>{</span>
    <span>temp_mat</span><span>.</span><span>col</span><span>(</span><span>i</span><span>)</span> <span>=</span> <span>(</span><span>i</span> <span>%</span> <span>2</span> <span>?</span> <span>o_vec</span> <span>:</span> <span>e_vec</span><span>);</span> <span>// slight change
</span>  <span>}</span>
  
  <span>return</span> <span>temp_mat</span><span>;</span>
<span>}</span></code></pre>
</figure>
<h4>Calculation Engine</h4>
<p>Next, we need to create a computational loop to subset the appropriate continuous areas of the matrix to figure out the amount of neighbors. In comparison to the problem statement, note that this loop is without the <code>t</code> as we no longer need to repeat calculations within this approach. Instead, we only need to compute the values once and then cache the result before duplicating it across the array.</p>
<figure>
<pre><code data-lang="cpp"><span>#include &lt;RcppArmadillo.h>
</span><span>// [[Rcpp::depends(RcppArmadillo)]]
</span>
<span>// --- Calculation engine
</span>
<span>// [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>mat</span> <span>calc_matrix</span><span>(</span><span>arma</span><span>::</span><span>mat</span> <span>temp_mat</span><span>){</span>
  
  <span>unsigned</span> <span>int</span> <span>xdim</span> <span>=</span> <span>temp_mat</span><span>.</span><span>n_rows</span><span>;</span>
  
  <span>unsigned</span> <span>int</span> <span>ydim</span> <span>=</span> <span>temp_mat</span><span>.</span><span>n_cols</span><span>;</span>
  
  <span>arma</span><span>::</span><span>mat</span> <span>res</span> <span>=</span> <span>temp_mat</span><span>;</span>
  
  <span>// Subset the rows
</span>  <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>x</span> <span>=</span> <span>2</span><span>;</span> <span>x</span> <span>&lt;</span> <span>xdim</span><span>-</span><span>2</span><span>;</span> <span>x</span><span>++</span><span>){</span>
    
    <span>arma</span><span>::</span><span>mat</span> <span>temp_row_sub</span> <span>=</span> <span>temp_mat</span><span>.</span><span>rows</span><span>(</span><span>x</span><span>-</span><span>2</span><span>,</span> <span>x</span><span>+</span><span>2</span><span>);</span>
    
    <span>// Iterate over the columns with unit accumulative sum
</span>    <span>for</span> <span>(</span><span>unsigned</span> <span>int</span> <span>y</span> <span>=</span> <span>2</span><span>;</span> <span>y</span> <span>&lt;</span>  <span>ydim</span><span>-</span><span>2</span><span>;</span> <span>y</span><span>++</span><span>){</span>
      <span>res</span><span>(</span><span>x</span><span>,</span><span>y</span><span>)</span> <span>=</span> <span>accu</span><span>(</span><span>temp_row_sub</span><span>.</span><span>cols</span><span>(</span><span>y</span><span>-</span><span>2</span><span>,</span><span>y</span><span>+</span><span>2</span><span>));</span>
    <span>}</span>
  <span>}</span>
  
  <span>return</span> <span>res</span><span>;</span>
<span>}</span></code></pre>
</figure>
<h4>Call Main Function</h4>
<p><em>Whew</em>, that was a lot of work. But, by approaching the problem this way, we have:</p>
<ol>
<li>Created reusable code snippets.</li>
<li>Decreased the size of the main call function.</li>
<li>Improved clarity of each operation.</li>
</ol>
<p>Now, the we are ready to write the glue that combines all the different components. As a result, we will obtain the desired neighbor information.</p>
<figure>
<pre><code data-lang="cpp"><span>// --- Requires the previous code snippets to be within one file inorder to run.
</span>
<span>// --- Main Engine
</span>
<span>// Create the desired cube information
// [[Rcpp::export]]
</span><span>arma</span><span>::</span><span>cube</span> <span>dim_to_cube</span><span>(</span><span>unsigned</span> <span>int</span> <span>xdim</span> <span>=</span> <span>4</span><span>,</span> <span>unsigned</span> <span>int</span> <span>ydim</span> <span>=</span> <span>4</span><span>,</span> <span>unsigned</span> <span>int</span> <span>tdim</span> <span>=</span> <span>3</span><span>)</span> <span>{</span>
  
  <span>// Initialize values in A
</span>  <span>arma</span><span>::</span><span>cube</span> <span>res</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>);</span>
  
  <span>if</span><span>(</span><span>xdim</span> <span>%</span> <span>2</span> <span>==</span> <span>0</span><span>){</span>
    <span>res</span><span>.</span><span>each_slice</span><span>()</span> <span>=</span> <span>calc_matrix</span><span>(</span><span>make_even_matrix_case1</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>));</span>
  <span>}</span><span>else</span><span>{</span>
    
    <span>if</span><span>(</span><span>ydim</span> <span>%</span> <span>2</span> <span>==</span> <span>0</span><span>){</span>
      
      <span>res</span><span>.</span><span>each_slice</span><span>()</span> <span>=</span> <span>calc_matrix</span><span>(</span><span>make_odd_matrix_case3</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>));</span>
      
    <span>}</span><span>else</span><span>{</span>
      
      <span>arma</span><span>::</span><span>mat</span> <span>first_odd_mat</span> <span>=</span> <span>calc_matrix</span><span>(</span><span>make_odd_matrix_case2</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>));</span>
      
      <span>arma</span><span>::</span><span>mat</span> <span>sec_odd_mat</span> <span>=</span> <span>calc_matrix</span><span>(</span><span>make_odd_matrix_case3</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>));</span>
      
      <span>for</span><span>(</span><span>unsigned</span> <span>int</span> <span>t</span> <span>=</span> <span>0</span><span>;</span> <span>t</span> <span>&lt;</span> <span>tdim</span><span>;</span> <span>t</span><span>++</span><span>){</span>
        <span>res</span><span>.</span><span>slice</span><span>(</span><span>t</span><span>)</span> <span>=</span> <span>(</span><span>t</span> <span>%</span> <span>2</span> <span>?</span> <span>first_odd_mat</span>  <span>:</span> <span>sec_odd_mat</span><span>);</span>
      <span>}</span>
      
    <span>}</span>
    
  <span>}</span>
  
  <span>return</span> <span>res</span><span>;</span>
<span>}</span></code></pre>
</figure>
<h4>Verification of Results</h4>
<p>To verify, let’s quickly create similar cases and test them against the original R function:</p>
<h5>Case 1:</h5>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>6</span><span>;</span> <span>ydim</span> <span>=</span> <span>5</span><span>;</span> <span>tdim</span> <span>=</span> <span>3</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span><span>dim</span><span>=</span><span>c</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>));</span> <span>res</span> <span>=</span> <span>a</span>

<span>all.equal</span><span>(</span><span>dim_to_cube</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>),</span> <span>cube_r</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>[</span><span>1</span><span>]</span> <span>TRUE</span></code></pre>
</figure>
<h5>Case 2:</h5>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>7</span><span>;</span> <span>ydim</span> <span>=</span> <span>6</span><span>;</span> <span>tdim</span> <span>=</span> <span>3</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span><span>dim</span><span>=</span><span>c</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>));</span> <span>res</span> <span>=</span> <span>a</span>

<span>all.equal</span><span>(</span><span>dim_to_cube</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>),</span> <span>cube_r</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>[</span><span>1</span><span>]</span> <span>TRUE</span></code></pre>
</figure>
<h5>Case 3:</h5>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>7</span><span>;</span> <span>ydim</span> <span>=</span> <span>7</span><span>;</span> <span>tdim</span> <span>=</span> <span>3</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span><span>dim</span><span>=</span><span>c</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>));</span> <span>res</span> <span>=</span> <span>a</span>

<span>all.equal</span><span>(</span><span>dim_to_cube</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>),</span> <span>cube_r</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<figure>
<pre><code data-lang="r"><span>[</span><span>1</span><span>]</span> <span>TRUE</span></code></pre>
</figure>
<h2>Closing Time – You don’t have to go home but you can’t stay here.</h2>
<p>With all of these different methods now thoroughly described, let’s do one last benchmark to figure out the best of the best.</p>
<figure>
<pre><code data-lang="r"><span>xdim</span> <span>=</span> <span>20</span>
<span>ydim</span> <span>=</span> <span>20</span>
<span>tdim</span> <span>=</span> <span>5</span>
<span>a</span> <span>=</span> <span>array</span><span>(</span><span>0</span><span>:</span><span>1</span><span>,</span><span>dim</span><span>=</span><span>c</span><span>(</span><span>xdim</span><span>,</span><span>ydim</span><span>,</span><span>tdim</span><span>))</span>
<span>res</span> <span>=</span> <span>a</span>


<span>microbenchmark</span><span>::</span><span>microbenchmark</span><span>(</span><span>r_orig</span>      <span>=</span> <span>cube_r</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>),</span>
                               <span>r_cached</span>    <span>=</span> <span>cube_r_cache</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>),</span>
                               <span>cpp_core1</span>   <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>1</span><span>),</span> 
                               <span>cpp_core2</span>   <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>2</span><span>),</span> 
                               <span>cpp_core3</span>   <span>=</span> <span>cube_cpp_parallel</span><span>(</span><span>a</span><span>,</span> <span>res</span><span>,</span> <span>3</span><span>),</span>
                               <span>cpp_pattern</span> <span>=</span> <span>dim_to_cube</span><span>(</span><span>xdim</span><span>,</span> <span>ydim</span><span>,</span> <span>tdim</span><span>))</span></code></pre>
</figure>
<p><strong>Output:</strong></p>
<figure>
<pre><code data-lang="r"><span>Unit</span><span>:</span> <span>microseconds</span>
            <span>expr</span>      <span>min</span>        <span>lq</span>       <span>mean</span>    <span>median</span>        <span>uq</span>      <span>max</span> <span>neval</span>
          <span>r_orig</span> <span>3545.746</span> <span>3912.1155</span> <span>4220.26266</span> <span>4062.5115</span> <span>4270.3515</span> <span>6132.152</span>   <span>100</span>
        <span>r_cached</span> <span>2750.397</span> <span>3003.5730</span> <span>3338.90562</span> <span>3132.2870</span> <span>3299.1740</span> <span>5115.018</span>   <span>100</span>
       <span>cpp_core1</span>  <span>170.006</span>  <span>181.0935</span>  <span>195.04401</span>  <span>187.6615</span>  <span>202.9760</span>  <span>302.570</span>   <span>100</span>
       <span>cpp_core2</span>  <span>112.524</span>  <span>118.8390</span>  <span>135.98392</span>  <span>126.6020</span>  <span>142.9040</span>  <span>322.466</span>   <span>100</span>
       <span>cpp_core3</span>   <span>81.936</span>   <span>87.9210</span>  <span>120.86310</span>  <span>105.0160</span>  <span>150.9475</span>  <span>261.475</span>   <span>100</span>
     <span>cpp_pattern</span>   <span>40.240</span>   <span>44.0575</span>   <span>51.85476</span>   <span>48.7645</span>   <span>59.2875</span>   <span>76.207</span>   <span>100</span></code></pre>
</figure>
<p>As can be seen, the customized approach based on the data’s pattern provided the fastest speed up. Meanwhile, the code port into C++ yielded much better results than the optimized R version. Lastly, the parallelization within R was simply too time consuming for a one-off computation.</p>
<h2>Misc</h2>
<p>The benchmarks listed above were taken using <a target="_blank" href="https://cran.r-project.org/web/packages/microbenchmark/"><code>microbenchmark</code> package</a>.</p>
<p>Hardware:</p>
<ul>
<li>iMac (27-inch, Late 2013) with:
<ul>
<li>3.2 GHz Intel Core i5</li>
<li>32 GB 1600 MHz DDR3</li>
<li>NVIDIA GeForce GT 755M 1024MB</li>
</ul>
</li>
</ul>
<p>Software:</p>
<ul>
<li>R version 3.3.0</li>
<li>RcppArmadillo version 0.7.100.3.1</li>
<li>Rcpp version 0.12.5</li>
<li>Compiler: clang-omp under c++98</li>
</ul>
<p><img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/rcpp-gallery/~4/Fk0vitQyuQI"></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/optimizing-code-vs-recognizing-patterns-with-3d-arrays/"></div></div><div><a data-url="http://www.r-bloggers.com/optimizing-code-vs-recognizing-patterns-with-3d-arrays/" data-counturl="http://www.r-bloggers.com/optimizing-code-vs-recognizing-patterns-with-3d-arrays/" data-count="vertical" data-text="Optimizing Code vs Recognizing Patterns with 3D Arrays" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://gallery.rcpp.org//articles/optimizing-code-vs-recognizing-patterns-with-3d-arrays/"> Rcpp Gallery</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://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="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/ClUNDnIc6Uk"></a></div>