Python vectorizing nested for loops
<p>I'd appreciate some help in finding and understanding a pythonic way to optimize the following array manipulations in nested for loops:</p>

<pre><code>def _func(a, b, radius):
    "Return 0 if a&gt;b, otherwise return 1"
    if distance.euclidean(a, b) &lt; radius:
        return 1
    else:
        return 0

def _make_mask(volume, roi, radius):
    mask = numpy.zeros(volume.shape)
    for x in range(volume.shape[0]):
        for y in range(volume.shape[1]):
            for z in range(volume.shape[2]):
                mask[x, y, z] = _func((x, y, z), roi, radius)
    return mask
</code></pre>

<p>Where <code>volume.shape</code> (182, 218, 200) and <code>roi.shape</code> (3,) are both <code>ndarray</code> types; and <code>radius</code> is an <code>int</code> </p>

<p><strong>Approach #1</strong></p>

<p>Here's a vectorized approach -</p>

<pre><code>m,n,r = volume.shape
x,y,z = np.mgrid[0:m,0:n,0:r]
X = x - roi[0]
Y = y - roi[1]
Z = z - roi[2]
mask = X**2 + Y**2 + Z**2 &lt; radius**2
</code></pre>

<p>Possible improvement : We can probably speedup the last step with <code>numexpr</code> module -</p>

<pre><code>import numexpr as ne

mask = ne.evaluate('X**2 + Y**2 + Z**2 &lt; radius**2')
</code></pre>

<p><strong>Approach #2</strong></p>

<p>We can also gradually build the three ranges corresponding to the shape parameters and perform the subtraction against the three elements of <code>roi</code> on the fly without actually creating the meshes as done earlier with <code>np.mgrid</code>. This would be benefited by the use of <a href="http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html"><code>broadcasting</code></a> for efficiency purposes. The implementation would look like this -</p>

<pre><code>m,n,r = volume.shape
vals = ((np.arange(m)-roi[0])**2)[:,None,None] + \
       ((np.arange(n)-roi[1])**2)[:,None] + ((np.arange(r)-roi[2])**2)
mask = vals &lt; radius**2
</code></pre>

<p>Simplified version : Thanks to @Bi Rico for suggesting an improvement here as we can use <code>np.ogrid</code> to perform those operations in a bit more concise manner, like so -</p>

<pre><code>m,n,r = volume.shape    
x,y,z = np.ogrid[0:m,0:n,0:r]-roi
mask = (x**2+y**2+z**2) &lt; radius**2
</code></pre>

<hr>

<p><strong>Runtime test</strong></p>

<p>Function definitions -</p>

<pre><code>def vectorized_app1(volume, roi, radius):
    m,n,r = volume.shape
    x,y,z = np.mgrid[0:m,0:n,0:r]
    X = x - roi[0]
    Y = y - roi[1]
    Z = z - roi[2]
    return X**2 + Y**2 + Z**2 &lt; radius**2

def vectorized_app1_improved(volume, roi, radius):
    m,n,r = volume.shape
    x,y,z = np.mgrid[0:m,0:n,0:r]
    X = x - roi[0]
    Y = y - roi[1]
    Z = z - roi[2]
    return ne.evaluate('X**2 + Y**2 + Z**2 &lt; radius**2')

def vectorized_app2(volume, roi, radius):
    m,n,r = volume.shape
    vals = ((np.arange(m)-roi[0])**2)[:,None,None] + \
           ((np.arange(n)-roi[1])**2)[:,None] + ((np.arange(r)-roi[2])**2)
    return vals &lt; radius**2

def vectorized_app2_simplified(volume, roi, radius):
    m,n,r = volume.shape    
    x,y,z = np.ogrid[0:m,0:n,0:r]-roi
    return (x**2+y**2+z**2) &lt; radius**2
</code></pre>

<p>Timings -</p>

<pre><code>In [106]: # Setup input arrays  
     ...: volume = np.random.rand(90,110,100) # Half of original input sizes 
     ...: roi = np.random.rand(3)
     ...: radius = 3.4
     ...: 

In [107]: %timeit _make_mask(volume, roi, radius)
1 loops, best of 3: 41.4 s per loop

In [108]: %timeit vectorized_app1(volume, roi, radius)
10 loops, best of 3: 62.3 ms per loop

In [109]: %timeit vectorized_app1_improved(volume, roi, radius)
10 loops, best of 3: 47 ms per loop

In [110]: %timeit vectorized_app2(volume, roi, radius)
100 loops, best of 3: 4.26 ms per loop

In [139]: %timeit vectorized_app2_simplified(volume, roi, radius)
100 loops, best of 3: 4.36 ms per loop
</code></pre>

<p>So, as always <code>broadcasting</code> showing its magic for a crazy almost <strong><code>10,000x</code></strong> speedup over the original code and more than <strong><code>10x</code></strong> better than creating meshes by using on-the-fly broadcasted operations!</p>

python numpy vectorization