Numpy vectorized zero-order interpolation
<p>I have an array</p>

<pre><code>p = [[0.9, 0.95, 0.99],
     [0.89, 0.94, 0.98],
     [0.9, 0.95, 0.99],
     [0.91, 0.96, 0.97],
    ]
</code></pre>

<p>and a uniform random number for each row</p>

<pre><code>r = [0.5, 
     0.9,
     0.3,
     0.99] 
</code></pre>

<p>I want to know the last column index where p is still smaller than r, i.e.</p>

<pre><code>c = [0, 1, 0, 3]
</code></pre>

<p>For the last case, no value is smaller. If one places a column of 1s at the end of p, this would be 3. A -1 is also acceptable for me.</p>

<p>Naive solution:</p>

<pre><code>c = []
for prow, ri in zip(p, r):
    ci = numpy.interp(ri, prow, arange(len(prow)))
    c.append(int(numpy.ceil(ci)))

print c
[0, 1, 0, 3]
</code></pre>

<p>But I am looking for a vectorised solution that is fast and works for large arrays (millions of rows, ~10 columns).</p>

<p>I looked into these solutions:</p>

<ul>
<li>scipy.interpolate.interp1d(kind=zero)

<ul>
<li>this seems to require a outer python loop</li>
</ul></li>
<li>r > p and numpy.where

<ul>
<li>also seems to require a outer python loop</li>
</ul></li>
<li>using numpy.random.choice</li>
</ul>

<p>For the last one I would place (differential) probabilities instead of cumulative ones:</p>

<pre><code>p = [[0.9, 0.05, 0.04],
     [0.89, 0.05, 0.04],
     [0.9, 0.05, 0.04],
     [0.91, 0.05, 0.01],
    ]
</code></pre>

<p>but numpy.random.choice does not support vectorization (<a href="https://github.com/numpy/numpy/issues/2724" rel="nofollow">1</a>, <a href="https://github.com/numpy/numpy/issues/5344" rel="nofollow">2</a>).</p>

<p>Is numpy.vectorise the solution, or Cython? I am looking for a fast solution.</p>

<p>Here's one vectorized solution using <a href="http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html" rel="nofollow"><code>broadcasting</code></a> -</p>

<pre><code>mask = (p &gt; r[:,None])
out = np.where(mask.any(1),mask.argmax(1),p.shape[1])
</code></pre>

<p>Sample run -</p>

<pre><code>In [50]: p
Out[50]: 
array([[ 0.9 ,  0.95,  0.99],
       [ 0.89,  0.94,  0.98],
       [ 0.9 ,  0.95,  0.99],
       [ 0.91,  0.96,  0.97]])

In [51]: r
Out[51]: array([ 0.5 ,  0.9 ,  0.3 ,  0.99])

In [52]: mask = (p &gt; r[:,None]) # 2D Mask of row-wise comparisons

In [53]: mask
Out[53]: 
array([[ True,  True,  True],
       [False,  True,  True],
       [ True,  True,  True],
       [False, False, False]], dtype=bool)

In [54]: np.where(mask.any(1),mask.argmax(1),p.shape[1])
Out[54]: array([0, 1, 0, 3])
</code></pre>

python numpy optimization vectorization