numpy korelacji krzyżowej - wektoryzacji

Question

Mam wiele korelacji wzajemnych do obliczenia i szukam najszybszy sposób to zrobić. Zakładam, że wektoryzacja problemu pomogłaby raczej niż w przypadku pętli

Mam tablicę 3D oznaczoną jako próbka elektrody x punktualny x (kształt: 64x256x913). Chcę obliczyć maksymalną korelację krzyżową punktów czasowych dla każdej pary elektrod dla każdej próby.

W szczególności: dla każdej próby chcę wziąć każdą kombinację par elektrod i obliczyć maksymalną wartość korelacji krzyżowej dla każdej pary. To spowoduje 4096 (64 * 64) maks. Wartości korelacji krzyżowej w jednym wierszu/wektorze. Zostanie to zrobione dla każdej próby, układając każdy z rzędów/wektorów jeden na drugim, tworząc ostateczną tablicę 2D o kształcie 913 * 4096 zawierającą maksymalne wartości korelacji krzyżowej

To dużo obliczeń, ale chcę spróbuj znaleźć najszybszą metodę, aby to zrobić. Wyśmiałem trochę proto-kodu używając list jako kontenerów, które mogą nieco lepiej wyjaśnić problem. Mogą tam być pewne błędy logiczne, ale tak czy inaczej kod nie działa na moim komputerze, ponieważ jest tak dużo, aby obliczyć, że pyton po prostu zawiesza się. Oto ona:

#allData is a 64x256x913 array 

all_experiment_trials = [] 
for trial in range(allData.shape[2]): 
    all_trial_electrodes = [] 
    for electrode in range(allData.shape[0]): 
     for other_electrode in range(allData.shape[0]): 
      if electrode == other_electrode: 
       pass 
      else: 
       single_xcorr = max(np.correlate(allData[electrode,:,trial], allData[other_electrode,:,trial], "full")) 
       all_trial_electrodes.append(single_xcorr) 
    all_experiment_trials.append(all_trial_electrodes) 

Oczywiście pętle są bardzo powolne dla tego typu rzeczy. Czy istnieje wektoryzowane rozwiązanie tego przy użyciu numpy array?

Mam wyrejestrowany rzeczy jak correlate2d() i tym podobne, ale nie sądzę, że naprawdę działa w moim przypadku tak nie jestem pomnożenie 2 matryce wraz

Answer 1

Oto jeden wektorowy podejście oparte na np.einsum -

def vectorized_approach(allData): 
    # Get shape 
    M,N,R = allData.shape 

    # Valid mask based on condition: "if electrode == other_electrode" 
    valid_mask = np.mod(np.arange(M*M),M+1)!=0 

    # Elementwise multiplications across all elements in axis=0, 
    # and then summation along axis=1 
    out = np.einsum('ijkl,ijkl->lij',allData[:,None,:,:],allData[None,:,:,:]) 

    # Use valid mask to skip columns and have the final output 
    return out.reshape(R,-1)[:,valid_mask] 

Runtime testów i weryfikacji wyników -

In [10]: allData = np.random.rand(20,80,200) 

In [11]: def org_approach(allData): 
    ...:  all_experiment_trials = [] 
    ...:  for trial in range(allData.shape[2]): 
    ...:   all_trial_electrodes = [] 
    ...:   for electrode in range(allData.shape[0]): 
    ...:    for other_electrode in range(allData.shape[0]): 
    ...:     if electrode == other_electrode: 
    ...:      pass 
    ...:     else: 
    ...:      single_xcorr = max(np.correlate(allData[electrode,:,trial], allData[other_electrode,:,trial])) 
    ...:      all_trial_electrodes.append(single_xcorr) 
    ...:   all_experiment_trials.append(all_trial_electrodes) 
    ...:  return all_experiment_trials 
    ...: 

In [12]: %timeit org_approach(allData) 
1 loops, best of 3: 1.04 s per loop 

In [13]: %timeit vectorized_approach(allData) 
100 loops, best of 3: 15.1 ms per loop 

In [14]: np.allclose(vectorized_approach(allData),np.asarray(org_approach(allData))) 
Out[14]: True