How to programm a stencil with Dask
<p>In many occasions, scientists simulates a system's dynamics using a Stencil, this is convolving a mathematical operator over a grid. Commonly, this operation consumes a lot of computational resources. <a href="https://en.wikipedia.org/wiki/Stencil_code" rel="nofollow">Here</a> is a good explanation of the idea. </p>

<p>In numpy, the canonical way of programming a 2D 5-points stencil is as follows:</p>

<pre><code>for i in range(rows):
    for j in range(cols):
        grid[i, j] = ( grid[i,j] + grid[i-1,j] + grid[i+1,j] + grid[i,j-1] + grid[i,j+1]) / 5
</code></pre>

<p>Or, more efficiently, using slicing: </p>

<pre><code>grid[1:-1,1:-1] = ( grid[1:-1,1:-1] + grid[0:-2,1:-1] + grid[2:,1:-1] + grid[1:-1,0:-2] + grid[1:-1,2:] ) / 5
</code></pre>

<p>However, if your grid is really big, it won't fix in your memory, or if the convolution operation is really complicated it will take a very long time, parallel programing techniques are use to overcome this problems or simply to get the result faster. Tools like <a href="http://dask.pydata.org/en/latest/" rel="nofollow">Dask</a> allow scientist to program this simulations by themselves, in a parallel-almost-transparent manner. Currently, Dask doesn't support item assignment, so, how can I program a stencil with Dask. </p>

<p>Dask internally divides arrays into smaller numpy arays, when you create an array with dask.array, you must provide some information about how to divide it into <em>chunks</em>, like this:</p>

<pre><code>grid = dask.array.zeros((100,100), chunks=(50,50))
</code></pre>

<p>That requests an array of 100 x 100 divided in 4 chunks. Now, to convolve an operation over the newly created array, information of chunk's borders must be shared. <a href="http://dask.pydata.org/en/latest/array-ghost.html" rel="nofollow">Dask ghost cells</a>, manages situations like this. </p>

<p>A common work-flow implies:</p>

<ol>
<li>Creating the array (if it didn't exist before)</li>
<li>Commanding the creation of ghost cells</li>
<li>Mapping a calculation</li>
<li>Trimming the borders</li>
</ol>

<p>For example,</p>

<pre><code>import dask.array as da
grid = da.zeros((100,100), chunks=(50,50))
g = da.ghost.ghost(grid, depth={0:1,1:1}, boundary={0:0,1:1})
g2 = g.map_blocks( some_function ) 
s = da.ghost.trim_internals(g2, {0:1,1:1})
s.compute()
</code></pre>

<p>Remember that Dask creates a dictionary to represent the task graph, the real computation is triggered by <code>s.compute()</code>. As noted by MRocklin, the mapped function must return a numpy array.  </p>

<h2>A note about schedulers</h2>

<p>By default, dask.array uses dask.theated scheduler to boost performance, but once the information is shared, problems similar to a stencil are embarrassing parallel, that means no resources nor information must be shared, and computations can be mapped to different cores even different computers. To do this, one can instruct dask to use a different scheduler, for example dask.multiprocessing:</p>

<pre><code>import dask.multiprocessing
import dask

dask.set_options(get=dask.multiprocessing.get)
</code></pre>

<p>When <code>compute()</code> is triggered, Dask will create multiple instances of python, if your applications is big enough to pay the overhead of creating this new instances, maybe dask.multiprocessing will deliver a better performance. 
More information about Dask schedulers can be found <a href="http://dask.pydata.org/en/latest/scheduler-overview.html" rel="nofollow">here</a>.</p>

python dask