Search code examples
pythonarraysnumpyrandomvectorization

generate a 2D array of numpy.random.choice without replacement


I'm tyring to make my code faster by removing some for loops and using arrays. The slowest step right now is the generation of the random lists.

context: I have a number of mutations in a chromosome, i want to perform 1000 random "chromosomes" with the same length and same number of mutation but their positions are randomized.

here is what I'm currently running to generate these randomized mutation positions:

iterations=1000
Chr_size=1000000
num_mut=500
randbps=[]
for k in range(iterations):
    listed=np.random.choice(range(Chr_size),num_mut,replace=False)
    randbps.append(listed)

I want to do something similar to what they cover in this question

np.random.choice(range(Chr_size),size=(num_mut,iterations),replace=False)

however without replacement applies to the array as a whole.

further context: later in the script i go through each randomized chromosome and count the number of mutations in a given window:

for l in range(len(randbps)):

    arr=np.asarray(randbps[l])

    for i in range(chr_last_window[f])[::step]:
    
        counter=((i < arr) & (arr < i+window)).sum()

Solution

  • I don't know how np.random.choice is implemented but I am guessing it is optimized for a general case. Your numbers, on the other hand, are not likely to produce the same sequences. Sets may be more efficient for this case, building from scratch:

    import random
    
    def gen_2d(iterations, Chr_size, num_mut):
        randbps = set()
        while len(randbps) < iterations:
            listed = set()
            while len(listed) < num_mut:
                listed.add(random.choice(range(Chr_size)))
            randbps.add(tuple(sorted(listed)))
        return np.array(list(randbps))
    

    This function starts with an empty set, generates a single number in range(Chr_size) and adds that number to the set. Because of the properties of the sets it cannot add the same number again. It does the same thing for the randbps as well so each element of randbps is also unique.

    For only one iteration of np.random.choice vs gen_2d:

    iterations=1000
    Chr_size=1000000
    num_mut=500
    
    %timeit np.random.choice(range(Chr_size),num_mut,replace=False)
    10 loops, best of 3: 141 ms per loop
    
    %timeit gen_2d(1, Chr_size, num_mut)
    1000 loops, best of 3: 647 µs per loop