# Fastest way for working with itertools.combinations

I need to speed up the function below:

import numpy as np
import itertools

def combcol(myarr):
ndims = myarr.shape[0]
solutions = []
for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(np.arange(ndims), 6):
c1, c2, c3, c4, c5, c6 = myarr[idx1,1], myarr[idx2,2], myarr[idx3,1], myarr[idx4,2], myarr[idx5,1], myarr[idx6,2]
if c1-c2>0 and c2-c3<0 and c3-c4>0 and c4-c5<0 and c5-c6>0 :
solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
return solutions

X = np.random.random((20, 10))
Y = np.random.random((40, 10))

if __name__=='__main__':
from timeit import Timer
t = Timer(lambda : combcol(X))
t1 = Timer(lambda : combcol(Y))
print('t : ',t.timeit(number=1),'t1 : ',t1.timeit(number=1))


Results:

t :  0.6165180211451455 t1 :  64.49216925614847


The algorithm is too slow for my standard use (myarr.shape[0] = 500). Is there a NumPy way to decrease execution time of this function (without wasting too much memory)? Is it possible to implement the problem in Cython?

-
What are you trying to achieve here? –  Gareth Rees Dec 29 '13 at 20:08

Your algorithm is of time complexity O(n^6) in ndims. If your current code needs 65 s for ndims = 40, then one could estimate the time for ndims = 500 to be roughly 8 years. Even an improvement of the most inner block's execution time by a factor of 1000 would yield an execution time of about 3 days. I guess this is not good enough.

So you should better try to improve the algorithm itself. However I do not really understand the problem the function is trying to solve and I couldn't come up with a better solution on the spot.

Also I haven't looked at numpy approaches, because you seem to have no vectorizable calculations in there.

Nonetheless I tried to make improvements to the code as is. I came up with the following:

1. You use only the a small fraction of myarr, so extract these parts first

myarr1 = myarr[:,1]
myarr2 = myarr[:,2]

2. Since you do many index lookups on these arrays turn them into lists. Numpy is not good at single index lookups.

myarr1 = myarr[:,1].tolist()
myarr2 = myarr[:,2].tolist()

3. Instead of substracting and comparing to 0, compare the variables directly saving one operation:

if c1>c2 and c2<c3 and c3>c4 and c4<c5 and c5>c6:

4. Instead of using and, you can also connect the comparisons all together, potentially saving the double variable lookups:

 if c1 > c2 < c3 > c4 < c5 > c6:

5. Instead of assigning the array elements to c1, c2, etc., use them directly. Since later comparisons are not executed if earlier one already return False this might save array lookups:

if myarr1[idx1] > myarr2[idx2] < myarr1[idx3] > myarr2[idx4] < myarr1[idx5] > myarr2[idx6]:


All in all I get the following code:

import numpy as np
import itertools

def combcol(myarr):
ndims = myarr.shape[0]
solutions = []
for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(np.arange(ndims), 6):
c1, c2, c3, c4, c5, c6 = myarr[idx1,1], myarr[idx2,2], myarr[idx3,1], myarr[idx4,2], myarr[idx5,1], myarr[idx6,2]
if c1-c2>0 and c2-c3<0 and c3-c4>0 and c4-c5<0 and c5-c6>0 :
solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
return solutions

def combcol2(myarr):
ndims = myarr.shape[0]
myarr1 = myarr[:,1].tolist()
myarr2 = myarr[:,2].tolist()
solutions = []
for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(range(ndims), 6):
if myarr1[idx1] > myarr2[idx2] < myarr1[idx3] > myarr2[idx4] < myarr1[idx5] > myarr2[idx6]:
solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(myarr1[idx1], myarr2[idx2], myarr1[idx3], myarr2[idx4], myarr1[idx5], myarr2[idx6])))
return solutions

X = np.random.random((40, 10))
Y = [np.random.random((10, 10)) for i in range(100)]

if __name__=='__main__':
from timeit import Timer
for test_case in Y: assert combcol(test_case) == combcol2(test_case)
t = Timer(lambda : combcol(X))
t1 = Timer(lambda : combcol2(X))
print('t : ',t.timeit(number=1),'t1 : ',t1.timeit(number=1))


With results like these:

t :  12.221544027328491 t1 :  0.8104081153869629
t :  10.69440507888794 t1 :  0.6592199802398682
t :  12.91178011894226 t1 :  0.8727250099182129
t :  11.21804690361023 t1 :  0.6851639747619629


The improvement is not even close to what you need. Also I tried to run the algorithm with ndims = 80 and got memory problems. The reason is the number of solutions: There were 18679356 solutions taking up 5 GB of RAM. With other values for ndims I found that around 5% - 25% of all possible solutions are valid! So it seems you need to change the way you output your solutions if you want to get beyond ndims = 100 on any realistic machine.

Edit: I just notice that 5 GB is a bit too much for the number of solutions. However this was the memory allocated by the python process. A run with ndims = 100 was unsuccessful due to lack of available memory.

Edit2: There is actually something far more efficient you can do, by interweaving your solution checks with the combination generation (abandoning itertools.combinations). This doesn't look so nice however:

import numpy as np
import itertools

def combcol(myarr):
ndims = myarr.shape[0]
solutions = []
for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(np.arange(ndims), 6):
c1, c2, c3, c4, c5, c6 = myarr[idx1,1], myarr[idx2,2], myarr[idx3,1], myarr[idx4,2], myarr[idx5,1], myarr[idx6,2]
if c1-c2>0 and c2-c3<0 and c3-c4>0 and c4-c5<0 and c5-c6>0 :
solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
return solutions

def combcol2(myarr):
ndims = myarr.shape[0]
myarr1 = myarr[:,1].tolist()
myarr2 = myarr[:,2].tolist()
solutions = []

for idx1 in range(ndims):
c1 = myarr1[idx1]
for idx2 in range(idx1+1, ndims):
c2 = myarr2[idx2]
if c1 > c2:
for idx3 in range(idx2+1, ndims):
c3 = myarr1[idx3]
if c2 < c3:
for idx4 in range(idx3+1, ndims):
c4 = myarr2[idx4]
if c3 > c4:
for idx5 in range(idx4+1, ndims):
c5 = myarr1[idx5]
if c4 < c5:
for idx6 in range(idx5+1, ndims):
c6 = myarr2[idx6]
if c5 > c6:
solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
return solutions

X = np.random.random((40, 10))
Y = [np.random.random((10, 10)) for i in range(100)]

if __name__=='__main__':
from timeit import Timer
for test_case in Y: assert combcol(test_case) == combcol2(test_case)
t = Timer(lambda : combcol(X))
t1 = Timer(lambda : combcol2(X))
print('t : ',t.timeit(number=1),'t1 : ',t1.timeit(number=1))


The result:

t :  11.706523180007935 t1 :  0.15908312797546387
t :  12.531341075897217 t1 :  0.16385602951049805
t :  11.581043004989624 t1 :  0.14388418197631836


So with that we already have nearly a factor of 100x speed-up.

-
Thank you very much. –  querzy Dec 30 '13 at 13:12