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>Nice question.  You're correct that <a href="http://dask.pydata.org/en/latest/array.html" rel="nofollow">dask.array</a> <em>do</em> provide parallel computing but <em>don't</em> doesn't support item assignment.  We can solve stencil computations by making a function to operate on a block of numpy data at a time and then by mapping that function across our array with slightly overlapping boundaries.</p>

<h3>Pure Functions</h3>

<p>You should make a function that takes a numpy array and returns a new numpy array with the stencil applied.  This should not modify the original array.</p>

<pre><code>def apply_stencil(x):
    out = np.empty_like(x)
    ...  # do arbitrary computations on out    
    return out
</code></pre>

<h3>Map a function with overlapping regions</h3>

<p>Dask arrays operate in parallel by breaking an array into disjoint chunks of smaller arrays.  Operations like stencil computations will require a bit of overlap between neighboring blocks.  Fortunately this can be handled with the <a href="http://dask.pydata.org/en/latest/array-ghost.html" rel="nofollow">dask.array.ghost</a> module, and the <a href="http://dask.pydata.org/en/latest/array-api.html#dask.array.Array.map_overlap" rel="nofollow">dask.array.map_overlap</a> method in particular.</p>

<p>Actually, the example in the <code>map_overlap</code> docstring is a 1d forward finite difference computation</p>

<pre><code>&gt;&gt;&gt; x = np.array([1, 1, 2, 3, 3, 3, 2, 1, 1])
&gt;&gt;&gt; x = from_array(x, chunks=5)
&gt;&gt;&gt; def derivative(x):
...     return x - np.roll(x, 1)
&gt;&gt;&gt; y = x.map_overlap(derivative, depth=1, boundary=0)
&gt;&gt;&gt; y.compute()
array([ 1,  0,  1,  1,  0,  0, -1, -1,  0])
</code></pre>

python dask