Оптимизация операции отображения массива в python - PullRequest
Новая
       36

Оптимизация операции отображения массива в python

1 голос
/ 08 апреля 2020

Я пытаюсь избавиться от неэффективного набора вложенных циклов в python. У меня есть массив, который я назову S (f k , f q ), который необходимо сопоставить с другим массивом, который я назову Z (f i , α J ). Аргументы все частоты дискретизации. Оба массива имеют одинаковые размеры, которые выбираются пользователем. Правило отображения довольно простое:

f i = 0,5 · (f k - f q )
α j = f k + f q

В настоящее время я выполняю это с помощью ряда вложенных циклов: 

import numpy as np
nrows = 64
ncolumns = 16384
fk = np.fft.fftfreq(nrows)
fq = np.fft.fftfreq(ncolumns)
# using random numbers here to simplify the example
# in practice S is the result of several FFTs and complex multiplications
S = np.random.random(size=(nrows,ncolumns)) + 1j*np.random.random(size=(nrows,ncolumns))

fi = []
alphaj = []
Z = []
for k in range(-nrows//2,nrows//2):
    for q in range(-ncolumns//2,ncolumns//2):
        fi.append(0.5*(fk[k] - fq[q]))
        alphaj.append(fk[k] + fq[q])
        Z.append(S[k,q])

Очевидно, что это крайне неэффективно - при таком подходе операция отображения занимает больше времени, чем фактическое вычисление S (что на практике является результатом нескольких БПФ и сложных умножений). Я хотел бы найти способ векторизовать это, но у меня возникают проблемы с правильным подходом. Будем весьма благодарны за любые предложения.

Примечание: это связано с другим вопросом о том, как сохранить результаты. Поскольку речь идет об оптимизации, я подумал, что было бы лучше создать два отдельных вопроса.

Ответы [ 2 ]

0 голосов
/ 08 апреля 2020

Это не использует отрицательное индексирование вашей исходной функции, но, возвращая массивы, вы можете использовать нормальное индексирование для отображения значений 

def weirdMath():
    nrows = 64
    ncolumns = 16384
    fk = np.fft.fftfreq(nrows)
    fq = np.fft.fftfreq(ncolumns)
    S = np.random.random(size=(nrows,ncolumns)) + 1j*np.random.random(size=(nrows,ncolumns))

    fi = .5*fk[:,np.newaxis] - fq
    alphaj = fk[:,np.newaxis] + fq
    return fi, alphaj, S


>>> f1,a1=weirdMath()
>>> f1.size
1048576
>>> f1[32,:10]
array([ 0.25      ,  0.24993896,  0.24987793,  0.24981689,  0.24975586,
        0.24969482,  0.24963379,  0.24957275,  0.24951172,  0.24945068])

Доказательство с добавлением вращения осей для соответствия порядку вывода в исходном коде , Примечание: S был изменен на np.arange (), чтобы сравнение значений между функциями могло быть напрямую сопоставлено: 

def origCode():
    nrows = 64
    ncolumns = 16384
    fk = np.fft.fftfreq(nrows)
    fq = np.fft.fftfreq(ncolumns)
# using random numbers here to simplify the example
# in practice S is the result of several FFTs and complex multiplications
    #S = np.random.random(size=(nrows,ncolumns)) + 1j*np.random.random(size=(nrows,ncolumns))
    S = np.arange(nrows*ncolumns).reshape(nrows, ncolumns)
    fi = []
    alphaj = []
    Z = []

    for k in range(-nrows//2,nrows//2):
        for q in range(-ncolumns//2,ncolumns//2):
            fi.append(0.5*fk[k] - fq[q])
            alphaj.append(fk[k] + fq[q])
            Z.append(S[k,q])
    return fi, alphaj,Z 


def weirdMathWithRoll():
    nrows = 64
    ncolumns = 16384
    rowRollAdj = nrows%2
    colRollAdj = ncolumns%2

    fk = np.roll(np.fft.fftfreq(nrows), shift=(-nrows//2) + rowRollAdj, axis=0)

    fq = np.roll(np.fft.fftfreq(ncolumns), (-ncolumns//2) + colRollAdj)
    S = np.random.random(size=(nrows,ncolumns)) + 1j*np.random.random(size=(nrows,ncolumns))
    S = np.arange(nrows*ncolumns).reshape(nrows, ncolumns)
    s2 = np.roll(S,ncolumns//2 + colRollAdj, axis=1)
    s3 = np.roll(s2,nrows//2 + rowRollAdj, axis=0)

    fi = .5*fk[:,np.newaxis] - fq
    alphaj = fk[:,np.newaxis] + fq

    return fi, alphaj, s3

def testMath():
    f,a,z = origCode()
    f1,a1,s1 = weirdMathWithRoll()

    fMatch = f==f1.flatten()
    aMatch = a==a1.flatten()
    sMatch = z==s1.flatten()
    return (np.all(fMatch), np.all(aMatch), np.all(sMatch))

Вывод доказательства: 

>>> testMath()
(True, True, True)

Улучшение производительности: 

>>> timeit.timeit(origCode, number = 1)
0.984715332997439
>>> timeit.timeit(weirdMathWithRoll, number = 1)
0.051891374998376705
0 голосов
/ 08 апреля 2020

Индексирование с отрицательными значениями k делает то, что вы хотите? В Python / numpy fk [-1] означает последний, fk [-2] означает второй после последнего, et c. 

In [90]: S = np.arange(1,11)                                                                           
In [91]: Z = []                                                                                        
In [92]: for k in range(-5,5): 
    ...:     Z.append(S[k]) 
    ...:                                                                                               

In [94]: S                                                                                             
Out[94]: array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
In [95]: Z                                                                                             
Out[95]: [6, 7, 8, 9, 10, 1, 2, 3, 4, 5]

Или с нарезкой: 

In [96]: np.concatenate([S[5:],S[:5]])                                                                 
Out[96]: array([ 6,  7,  8,  9, 10,  1,  2,  3,  4,  5])
...