Использование pyfftw должным образом для ускорения работы над numpy

Question

Я посреди попыток сделать прыжок с Matlab на numpy, но мне отчаянно нужна скорость в моих fft. Теперь я знаю о pyfftw, но я не знаю, что я использую его правильно. Мой подход будет что-то вроде

import numpy as np 
import pyfftw 
import timeit 

pyfftw.interfaces.cache.enable() 

def wrapper(func, *args): 
    def wrapped(): 
     return func(*args) 
    return wrapped 

def my_fft(v): 
    global a 
    global fft_object 
    a[:] = v 
    return fft_object() 

def init_cond(X): 
    return my_fft(2.*np.cosh(X)**(-2)) 

def init_cond_py(X): 
    return np.fft.fft(2.*np.cosh(X)**(-2)) 

K = 2**16 
Llx = 10. 
KT = 2*K 
dx = Llx/np.float64(K) 
X = np.arange(-Llx,Llx,dx) 

global a 
global b 
global fft_object 
a = pyfftw.n_byte_align_empty(KT, 16, 'complex128') 
b = pyfftw.n_byte_align_empty(KT, 16, 'complex128') 
fft_object = pyfftw.FFTW(a,b) 

wrapped = wrapper(init_cond, X) 
print min(timeit.repeat(wrapped,repeat=100,number=1)) 

wrapped_two = wrapper(init_cond_py, X) 
print min(timeit.repeat(wrapped_two,repeat=100,number=1)) 

Я понимаю, что есть модифицирующие функции, а также стандартные интерфейсы к SciPy и NumPy FFT звонки через pyfftw. Они все вели себя очень медленно, хотя. Сначала создав экземпляр fft_object, а затем используя его глобально, я смог получить скорость так же быстро или немного быстрее, чем вызов fty numpy.

Это, как говорится, я работаю в предположении, что мудрость неявно хранится. Это правда? Мне нужно сделать это явным? Если да, то каков наилучший способ сделать это?

Кроме того, я думаю, что timeit полностью непрозрачен. Правильно ли я использую его? Он хранит мудрость, как я повторяю? Заранее спасибо за любую помощь, которую вы могли бы дать.

Answer 1

В интерактивной (IPython) сессии, я думаю, что следующее, что вы хотите сделать (timeit очень хорошо обрабатывается IPython):

In [1]: import numpy as np 

In [2]: import pyfftw 

In [3]: K = 2**16 

In [4]: Llx = 10. 

In [5]: KT = 2*K 

In [6]: dx = Llx/np.float64(K) 

In [7]: X = np.arange(-Llx,Llx,dx) 

In [8]: a = pyfftw.n_byte_align_empty(KT, 16, 'complex128') 

In [9]: b = pyfftw.n_byte_align_empty(KT, 16, 'complex128') 

In [10]: fft_object = pyfftw.FFTW(a,b) 

In [11]: a[:] = 2.*np.cosh(X)**(-2) 

In [12]: timeit np.fft.fft(a) 
100 loops, best of 3: 4.96 ms per loop 

In [13]: timeit fft_object(a) 
100 loops, best of 3: 1.56 ms per loop 

In [14]: np.allclose(fft_object(a), np.fft.fft(a)) 
Out[14]: True 

Вы читали tutorial? Что ты не понимаешь?

Я бы рекомендовал использовать объект builders interface для создания объекта FFTW. Имейте игру с различными настройками, самое главное количество потоков.

Мудрость не сохраняется по умолчанию. Вам нужно extract it yourself.

Все ваши globals не нужны - объекты, которые вы хотите изменить, изменяемы, поэтому вы можете справиться с ними просто отлично. fft_object всегда указывает на одно и то же, поэтому никаких проблем с тем, чтобы не быть глобальным. В идеале вы просто не хотите, что петля над ii. Я предлагаю, как структурировать ваши массивы, чтобы вы могли делать все ваши операции за один раз.

Редактировать: Редактировать: Я написал следующий абзац, только с беглым взглядом на ваш код и четко с это рекурсивное обновление, векторизация - не очевидный подход без какой-либо серьезной хитрости. У меня есть несколько комментариев о вашей реализации внизу.] Я подозреваю, что ваша проблема - более фундаментальное непонимание того, как наилучшим образом использовать язык, например Python (или даже Matlab) для цифровой обработки. Основной принцип как можно более векторный. Под этим я подразумеваю, что ваши вызовы python должны быть как можно меньше. Я не вижу, как это сделать с вашим примером, к сожалению (хотя я думал об этом только 2 минуты). Если это все еще не удается, подумайте о cython - хотя убедитесь, что вы действительно хотите спуститься по этому маршруту (т. Е. Вы исчерпали другие варианты).

Что касается глобальных переменных: не делайте этого. Если вы хотите создать объект с состоянием, используйте класс (для чего он предназначен) или, возможно, закрытие в вашем случае. Глобальное почти никогда не то, что вы хотите (я думаю, что у меня есть хотя бы смутно законное использование для него во всех моих письмах на python, и это в коде кеша в pyfftw). Предлагаю прочитать this nice SO question. Matlab - это дерьмовый язык - одна из многих причин для этого - его инструменты, которые, как правило, приводят к вредным привычкам.

Вам нужно всего лишь глобальное, если вы хотите изменить ссылку по всему миру.Я предлагаю прочитать немного больше о Python scoping rules и какие переменные really are в python.

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

Нижняя линия, не беспокойтесь об обновлении массивов при каждом вызове. Это почти наверняка не ваше узкое место, хотя убедитесь, что вы знаете о нормализации и отключите его, если хотите (это может немного замедлить ситуацию по сравнению с необработанным доступом методов update_arrays() и execute()).

Ваш код не использует кеш. Кэш используется только тогда, когда вы используете код interfaces и уменьшает накладные расходы Python при создании новых объектов FFTW внутри. Поскольку вы сами обрабатываете объект FFTW, нет никакой причины для кеша.

Код builders является менее ограниченным интерфейсом для получения объекта FFTW. Я почти всегда использую строителей сейчас (гораздо удобнее создавать объект FFTW с нуля). Случаи, в которых вы хотите создать объект FFTW напрямую, довольно редки, и мне было бы интересно узнать, что они собой представляют.

Комментарии к реализации алгоритма: Я не знаком с алгоритмом, который вы реализуете. Тем не менее, у меня есть несколько комментариев о том, как вы написали это на данный момент. Вы вычисляете nl_eval(wp) на каждом цикле, но насколько я могу сказать, это то же самое, что и nl_eval(w) из предыдущего цикла, поэтому вам не нужно его вычислять дважды (но это связано с предостережением, что довольно сложно посмотрите, что происходит, когда у вас есть глобальные перемены повсюду, поэтому я мог бы что-то упустить).

Не копируйте копии в my_fft или my_ifft. Просто сделайте fft_object(u) (2,29 мс против 1,67 мс на моей машине для переднего футляра). Процедура обновления внутреннего массива делает копию ненужной. Кроме того, как вы его написали, вы копируете дважды: c[:] означает «копировать в массив c», а массив, который вы копируете в c, составляет v.copy(), то есть копию v (всего две копии).

Более разумным (и, возможно, необходимым) является копирование вывода в массивы хранения (поскольку это позволяет избежать промежуточных результатов при вызове объекта FFTW), хотя убедитесь, что ваши массивы хранения правильно выровнены. Я уверен, что вы отметили, что это важно, но копировать результат довольно понятно.

Вы можете перемещать все ваши масштабирования вместе. 3 в вычислении wn можно перемещать внутри my_fft в nl_eval. Вы также можете объединить это с константой нормализации из ifft (и отключить его в pyfftw).

Посмотрите на numexpr для основных операций с массивами. Он может предложить немного ускорения по сравнению с ванильным количеством.

В любом случае возьмите то, что вы от всего этого.Без сомнения, я что-то пропустил или сказал что-то неправильное, поэтому, пожалуйста, примите его с таким смирением, какое я могу предложить. Стоит потратить немного времени на то, как Python клеится по сравнению с Matlab (на самом деле, просто забудьте о последнем).