Efficient use of numpy.random.choice with repeated numbers and alternatives
<p>I need to generate a large array with repeated elements, and my code is:</p>

<pre><code>np.repeat(xrange(x,y), data)
</code></pre>

<p>However, data is a numpy array with type float64 (but it represent integeres, no 2.1 there) and I get the error</p>

<pre><code>TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'
</code></pre>

<p>Exemple:</p>

<pre><code>In [35]: x
Out[35]: 26

In [36]: y
Out[36]: 50

In [37]: data
Out[37]: 
array([ 3269.,   106.,  5533.,   317.,  1512.,   208.,   502.,   919.,
     406.,   421.,  1690.,  2236.,   705.,   505.,   230.,   213.,
     307.,  1628.,  4389.,  1491.,   355.,   103.,   854.,   424.])
In [38]: np.repeat(xrange(x,y), data)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call    last)
&lt;ipython-input-38-105860821359&gt; in &lt;module&gt;()
----&gt; 1 np.repeat(xrange(x,y), data)

/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy    /core/fromnumeric.pyc in repeat(a, repeats, axis)
394         repeat = a.repeat
395     except AttributeError:
--&gt; 396         return _wrapit(a, 'repeat', repeats, axis)
397     return repeat(repeats, axis)
398 

/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy  /core/fromnumeric.pyc in _wrapit(obj, method, *args, **kwds)
 46     except AttributeError:
 47         wrap = None
---&gt; 48     result = getattr(asarray(obj), method)(*args, **kwds)
 49     if wrap:
 50         if not isinstance(result, mu.ndarray):

TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'
</code></pre>

<p>I solve it by changing the code to</p>

<pre><code>np.repeat(xrange(x,y), data.astype('int64'))
</code></pre>

<p>However, this is now one of the most expensive lines in my code!! Is there another alternative?</p>

<p>By the way, I using this inside</p>

<pre><code>np.random.choice(np.repeat(xrange(x,y), data.astype('int64')), z)
</code></pre>

<p>in order to get a sample without replacement with size z of the integers between x and y, with the number of each given in data. I guess this is the best approach for that also right?</p>

<h2>Problem statement</h2>

<p>Pretty interesting problem this one! Just to give the readers an idea about the problem without going into the minor data conversion issues, we have a range of values, let's say <code>a = np.arange(5)</code>, i.e.</p>

<pre><code>a = np.array([0,1,2,3,4])
</code></pre>

<p>Now, let's say we have another array with the number of repetitions listed for each of the <code>5</code> numbers in <code>a</code>. So, let those be :</p>

<pre><code>reps = np.array([2,4,6,2,2])
</code></pre>

<p>Next up, we are performing those repetitions :</p>

<pre><code>In [32]: rep_nums = np.repeat(a,reps)

In [33]: rep_nums
Out[33]: array([0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4])
</code></pre>

<p>Finally, we are looking to choose <code>z</code> number of elements out of those repeated numbers using <code>np.random.choice()</code> and without replacement.</p>

<p>Let's say <code>z = 7</code> to choose <code>7</code> elements, so with <code>np.random.choice()</code>, we would have :</p>

<pre><code>In [34]: np.random.choice(rep_nums,7,replace=False)
Out[34]: array([2, 4, 0, 2, 4, 1, 2])
</code></pre>

<p>Now, this <code>without replacement</code> term here might sound confusing, as we already have repeated numbers in <code>rep_nums</code>. But, what it essentially means is that, the output from <code>np.random.choice()</code> mustn't contain e.g. more than two <code>4's</code>, because <code>rep_nums</code> has two <code>4's</code>.</p>

<p>So, the problem is we want to get rid of that <code>np.repeat</code> part, which might be a bottleneck for really huge arrays.</p>

<h2>Proposed approach</h2>

<p>Looking at the output from <code>rep_nums</code>, one idea would be to generate <code>z = 7</code> unique elements ranging across the length of <code>rep_nums</code> :</p>

<pre><code>In [44]: np.random.choice(rep_nums.size,7,replace=False)
Out[44]: array([ 7,  2,  4, 10, 13,  8,  3])
</code></pre>

<p>These numbers represent indices for that length. So, we just need to look for the bin (out of the <code>5</code> bins) in <code>rep_nums</code> in which each of those <code>7</code>  numbers would go in. For that, we can use <code>np.searchsorted</code>. Thus, we would have an implementation to handle generic <code>x</code>, <code>y</code>, like so -</p>

<pre><code># Get the intervals of those bins
intervals = data.astype(int).cumsum()

# Decide length of array if we had repeated with `np.repeat`
max_num = intervals[-1]

# Get unique numbers (indices in this case)
ids = np.random.choice(max_num,z,replace=False)

# Use searchsorted to get bin IDs and add in `x` offset
out = x+np.searchsorted(intervals,ids,'right')
</code></pre>

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

<p>Functions :</p>

<pre><code>def org_app(x,y,z,data):
    rep_nums = np.repeat(range(x,y), data.astype('int64'))
    out = np.random.choice(rep_nums, z,replace=False)
    return out

def optimized_v1(x,y,z,data):     
    intervals = data.astype(int).cumsum()
    max_num = intervals[-1]
    ids = np.random.choice(max_num,z,replace=False)
    out = x+np.searchsorted(intervals,ids,'right')
    return out
</code></pre>

<p>Timings on full functions -</p>

<pre><code>In [79]: # Setup inputs
    ...: x = 100
    ...: y = 10010
    ...: z = 1000
    ...: data = np.random.randint(100,5000,(y-x)).astype(float)
    ...: 

In [80]: %timeit org_app(x,y,z,data)
1 loop, best of 3: 7.17 s per loop

In [81]: %timeit optimized_v1(x,y,z,data)
1 loop, best of 3: 6.92 s per loop
</code></pre>

<p>Doesn't look like we are getting good speedup. Let's dig deeper and find out how much are we saving on replacing <code>np.repeat</code>!</p>

<p>First off the original approach -</p>

<pre><code>In [82]: %timeit np.repeat(range(x,y), data.astype('int64'))
1 loop, best of 3: 227 ms per loop
</code></pre>

<p>Let's see how much improvement we got on this with the proposed approach. So, let's time everything except <code>np.random.choice()</code> in the proposed approach -</p>

<pre><code>In [83]: intervals = data.astype(int).cumsum()
    ...: max_num = intervals[-1]
    ...: ids = np.random.choice(max_num,z,replace=False)
    ...: out = x+np.searchsorted(intervals,ids,'right')
    ...: 

In [84]: %timeit data.astype(int).cumsum()
10000 loops, best of 3: 36.6 Âµs per loop

In [85]: %timeit intervals[-1]
10000000 loops, best of 3: 142 ns per loop

In [86]: %timeit x+np.searchsorted(intervals,ids,'right')
10000 loops, best of 3: 127 Âµs per loop
</code></pre>

<p><strong>This is much better than <code>227ms</code> from <code>np.repeat</code>!!</strong></p>

<p>So, we are hoping that at really huge arrays, the benefit from removing <code>np.repeat</code> would really shine, as otherwise <code>np.random.choice()</code> itself looks like the bottleneck.</p>

python python-2.7 numpy casting repeat