AGreatPigeon AGreatPigeon - 1 month ago 20
Python Question

Vectorizing numpy Multiple Condition Nested Loops

On attempting to produce Automatic Peak Detection in Noisy Periodic and Quasi-Periodic Signals, by Felix Scholkmann, Jens Boss and Martin Wolf in Python, I've hit a stumbling block in the implementation.

Upon attempting to optimise, I've noticed that the nested for loops are creating a bottleneck in processing time (taking 115394 ms on average to complete).
Is there a more efficient means of constructing the nested for loop?

N.B:
The parameter, signal, is a list of co-ordinates to which the algorithm will process which is of the form


-48701.0
-20914.0
-1757.0
-49278.0
-106781.0
-88139.0
-13587.0
28071.0
11880.0
-13375.0
-18056.0
-15248.0
-12476.0
-9832.0
-26365.0
-65734.0
-81657.0
-41566.0
6382.0
872.0
-30666.0
-20261.0
17543.0
6278.0
...


The list is 32768 lines long.

The function returns the indexes of the peaks detected to which is processed in another function.

def ampd(signal):

s_time = range(1, len(signal)+1)

[fitPolynomial, fitError] = np.polyfit(s_time, signal, 1)
fitSignal = np.polyval([fitPolynomial, fitError], s_time)

dtrSignal = signal - fitSignal

N = len(dtrSignal)
L = math.ceil(N/2.0)-1

creation_start = time.time()
np.random.seed(1969)

LSM = np.random.uniform(0, 2, size=(L, N))
creation_elapsedTime = time.time() - creation_start
print('LSM created in %s ms' % int(creation_elapsedTime * 1000))

loop_start = time.time()
for k in range(1, L):
for i in range(k+2, N-k+1):
if signal[i-1]>signal[i-k-1] and signal[i-1]>signal[i+k-1]:
LSM[k,i] = 0

loop_elapsedTime = time.time() - loop_start
print('Loop completed in %s ms' % int(loop_elapsedTime * 1000))

G = np.sum(LSM, axis=1)

l = min(enumerate(G), key=itemgetter(1))[0]

MLSM = LSM[0:l]

S = np.std(MLSM, ddof=1)

found_indices = np.where(MLSM == ((S-1) == 0))
del LSM
del MLSM

return found_indices[1]

JPG JPG
Answer

Here is a solution which uses only one loop

for k in range(1, L):
    mat=1-((signal[k+1:N-k]>signal[1:N-2*k]) & (signal[k+1:N-k]>signal[2*k+1:N]))
    LSM[k,k+2:N-k+1]*=mat

it's faster and seems do give the same solutions. You compare slices (as suggested by Ami Tavory), combine the comparisons with a &, which gives a True/False array; with 1-operation, you transform it to zeros and ones, the zeros corresponding to where the conditions are met. And lastly you multiply the row by the result.