Fast lookup from floating point index into array
<p>Suppose I have <code>n</code> intervals of length <code>len</code> on an axis. Each of them is assigned a value in an <code>ndarray</code>. Now I'd like to look-up some values at positions <code>query_pos</code> given as floating point numbers. Currently, I plan to do it the following way:</p>

<pre><code>n = 100
len = 0.39483
data = ... # creates ndarray with data of length n = 100
query_pos = ... # creates an ndarray with positions to query
values_at_query_pos = data[(query_pos / len).astype(int)] # questionable line
</code></pre>

<p><strong>Is that a good way to do it or are there more efficient ways to convert a floating point query position into an integer index and then read from an array?</strong> I particularly wonder, wether <code>astype(int)</code> is a cheap or an expensive operation, compared for example to the division or the memory read.</p>

<p>Some more remarks:</p>

<ul>
<li><p>Finally, it will be used in 2 and 3 dimensions. Currently, I plan to
catch positions that would lead to illegal indices <em>before</em> they go
into the lookup stage.</p></li>
<li><p>The <code>data</code> array will have a high-enough resolution so that I don't
need any filtering or interpolation. That will be done in previous
stages.</p></li>
</ul>

<p>Instead of dividing each element of <code>query_pos</code> by that scalar, we can pre-calculate the reciprocal of the scalar and use multiplication instead for some speedup there. The intuition is that division is a more costly affair than multiplication.</p>

<p>Here's a quick runtime test on it -</p>

<pre><code>In [177]: # Setup inputs
     ...: n = 100
     ...: len1 = 0.39483
     ...: data = np.random.rand(100)
     ...: query_pos = np.random.randint(0,25,(100000))
     ...: 

In [178]: %timeit query_pos / len1
1000 loops, best of 3: 612 Âµs per loop

In [179]: %timeit query_pos * (1/len1)
10000 loops, best of 3: 173 Âµs per loop
</code></pre>

<p>Secondly, if there are many repeated indices, just like in the setup used for the runtime test shown earlier, we can use <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.take.html" rel="nofollow"><code>np.take</code></a> for some further marginal improvement, as shown below -</p>

<pre><code>In [196]: %timeit data[(query_pos *(1/len1)).astype(int)]
1000 loops, best of 3: 538 Âµs per loop

In [197]: %timeit np.take(data,(query_pos * (1/len1)).astype(int))
1000 loops, best of 3: 526 Âµs per loop
</code></pre>

<p>If you are planning to use it on generic ndarrays, we would need to use <code>axis</code> param with <code>np.take</code>.</p>

<p>Comparing it with the original approach -</p>

<pre><code>In [202]: %timeit data[(query_pos / len1).astype(int)]
1000 loops, best of 3: 967 Âµs per loop
</code></pre>

<p>Finally, on the question of how the division operation stacks up against converting to <code>int</code>, they seem comparable on the big dataset. But, it seems you can't avoid the conversion as needed for indexing. Here's a timing test on it -</p>

<pre><code>In [210]: idx = query_pos * (1/len1)

In [211]: %timeit query_pos * (1/len1)
10000 loops, best of 3: 165 Âµs per loop

In [212]: %timeit idx.astype(int)
10000 loops, best of 3: 110 Âµs per loop
</code></pre>

python numpy optimization floating-point indices