numpy indexing with multiple arrays

Question

Given two sequences of data (of equal length) and quality values for each data point, I want to calculate a similarity score based upon a given scoring matrix.

What is the most efficient way to vectorize the following loop:

score = 0
for i in xrange(len(seq1)):
    score += similarity[seq1[i], seq2[i], qual1[i], qual2[i]]

similarity is a 4-dimensional float array, shape=(32, 32, 100, 100); seq1, seq2, qual1 and qual2 are 1-dimensional int arrays of equal length (of the order 1000 - 40000).

DSM · Accepted Answer · 2013-05-02 15:33:10Z

3

Shouldn't this Just Work(tm)?

>>> score = 0
>>> for i in xrange(len(seq1)):
        score += similarity[seq1[i], seq2[i], qual1[i], qual2[i]]
...     
>>> score
498.71792400493433
>>> similarity[seq1,seq2, qual1, qual2].sum()
498.71792400493433

Code:

import numpy as np

similarity = np.random.random((32, 32, 100, 100))
n = 1000
seq1, seq2, qual1, qual2 = [np.random.randint(0, s, n) for s in similarity.shape]

def slow():
    score = 0
    for i in xrange(len(seq1)):
        score += similarity[seq1[i], seq2[i], qual1[i], qual2[i]]
    return score

def fast():
    return similarity[seq1, seq2, qual1, qual2].sum()

gives:

>>> timeit slow()
100 loops, best of 3: 3.59 ms per loop
>>> timeit fast()
10000 loops, best of 3: 143 us per loop
>>> np.allclose(slow(),fast())
True

edited May 2, 2013 at 15:33

answered May 2, 2013 at 15:26

DSM

355k67 gold badges606 silver badges504 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

John Zwinck Over a year ago

Well that's just great. I like this better than my own answer (though I may leave mine just for contrast). +1.

GWW Over a year ago

That's a really neat numpy feature I didn't know about.

Quant Metropolis Over a year ago

Thanks! - Timings on my machine (including John Zwinck's answer as 'third') for 1000 iterations of length 10000: slow: 17.1289219856, fast: 0.61208987236, third: 15.7027080059

John Zwinck · Accepted Answer · 2013-05-02 15:31:13Z

0

How about this?

score = numpy.sum(map(similarity.__getitem__, zip(seq1, seq2, qual1, qual2)))

Of course you can try with itertools imap and izip too. The zip is necessary because __getitem__ takes a single tuple rather than four numbers...maybe that can be improved somehow by looking in a darker corner of the itertools module.

answered May 2, 2013 at 15:31

John Zwinck

252k44 gold badges346 silver badges459 bronze badges

Collectives™ on Stack Overflow

numpy indexing with multiple arrays

2 Answers 2

3 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related