Как размеры влияют на производительность в pyfftw? - PullRequest
Новая
       43

Как размеры влияют на производительность в pyfftw?

1 голос
/ 11 апреля 2019

Я пытаюсь реализовать трехмерную свертку, используя FFT с pyfftw. Я использовал в качестве базы код, размещенный в другом посте в SO: 

class CustomFFTConvolution(object):

def __init__(self, A, B, threads=1):

    shape = (np.array(A.shape) + np.array(B.shape))-1
    #shape=np.array(A.shape) - np.array(B.shape)+1
    if np.iscomplexobj(A) and np.iscomplexobj(B):
        self.fft_A_obj = pyfftw.builders.fftn(
                A, s=shape, threads=threads)
        self.fft_B_obj = pyfftw.builders.fftn(
                B, s=shape, threads=threads)
        self.ifft_obj = pyfftw.builders.ifftn(
                self.fft_A_obj.get_output_array(), s=shape,
                threads=threads)

    else:
        self.fft_A_obj = pyfftw.builders.rfftn(
                A, s=shape, threads=threads)
        self.fft_B_obj = pyfftw.builders.rfftn(
                B, s=shape, threads=threads)
        self.ifft_obj = pyfftw.builders.irfftn(
                self.fft_A_obj.get_output_array(), s=shape,
                threads=threads)

def __call__(self, A, B):
    s1=np.array(A.shape)
    s2=np.array(B.shape)

    fft_padded_A = self.fft_A_obj(A)
    fft_padded_B = self.fft_B_obj(B)

    ret= self.ifft_obj(fft_padded_A * fft_padded_B)

    return self._centered(ret, s1 - s2 + 1)

def _centered(self,arr, newshape):
    # Return the center newshape portion of the array.
    newshape = np.asarray(newshape)
    currshape = np.array(arr.shape)
    startind = (currshape - newshape) // 2
    endind = startind + newshape
    myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
    return arr[tuple(myslice)]

Мои данные A имеют форму (931, 411, 806), а мой фильтр B имеет форму (32, 32, 32). Если я запускаю этот код, используя 24 потока на машине с 24 ядрами, операция занимает 263 секунды. Теперь, если я проведу тот же эксперимент на той же машине, но на этот раз A имеет форму (806, 411, 931) , то есть просто изменение оси , код занимает всего 16 секунд. Что является причиной этого? Есть ли эмпирическое правило для получения наилучшей производительности? может быть, отступы одного из измерений? Спасибо!

1 Ответ

1 голос
/ 19 апреля 2019

Поскольку рассматривается дополнение, можно ли увеличить размер дополнения до четного или кратного небольшого числа простых чисел? Выбор четных размеров может разделить время настенных часов на 3.

В зависимости от размеров, некоторые алгоритмы DFT могут быть недоступны или неэффективны.Например, одним из наиболее эффективных алгоритмов выполнения ДПФ является алгоритм Кули-Тьюки .Он состоит в делении ДПФ сигнала составного размера N = N1 * N2 на N1 DTF размера N2.Как следствие, работает лучше для составных размеров, полученных умножением малых простых множителей (2, 3, 5, 7) , для которых в FFTW предусмотрены специальные эффективные алгоритмы.Из документации FFTW :

Например, стандартный дистрибутив FFTW работает наиболее эффективно для массивов, размер которых можно разбить на небольшие простые числа (2, 3, 5 и 7).), а в остальном он использует более медленную процедуру общего назначения.Если вам нужны эффективные преобразования других размеров, вы можете использовать генератор кода FFTW, который создает быстрые C-программы («кодлеты») для любого конкретного размера массива, который вас может заинтересовать.Например, если вам нужны преобразования размером 513 = 19 * 33, вы можете настроить FFTW для эффективной поддержки фактора 19.

Ваши размеры дополнения имеют высокие простые коэффициенты: 

931=>962=2*13*37
411=>442=2*13*17
806=>837=3*3*3*31

Заполнение может быть расширено, чтобы приблизиться к числам с небольшими простыми числами, например, например, 980, 448 и 864,Тем не менее, заполнение трехмерного изображения приводит к значительному увеличению объема памяти, до такой степени, что это не всегда возможно.

Почему изменение порядка размеров меняет время вычисления?Разница может быть связана с тем, что входной массив является действительным. Следовательно, DFT R2C выполняется над одним из измерений, затем C2C над вторым и третьим для вычисления 3D DFT.Если размер первого измерения, подлежащего преобразованию, является четным, преобразование R2C можно превратить в сложное ДПФ, равное половине размера, как показано здесь .Этот трюк не работает для нечетного размера.Как следствие, некоторые быстрые алгоритмы, вероятно, станут доступны после переворота 962 и 837.

Вот код для проверки: 

import pyfftw
import matplotlib.pyplot as plt
import multiprocessing
import numpy as np
from timeit import default_timer as timer

def listofgoodsizes():
    listt=[]
    p2=2
    for i2 in range(11):
        p3=1
        for i3 in range(7):
            p5=1
            for i5 in range(2):

                listt.append(p2*p3*p5)
                p5*=5
            p7=1
            for i7 in range(2):
                listt.append(p2*p3*p7)
                p7*=7

            p3*=3
        p2*=2
    listt.sort()
    return listt

def getgoodfftwsize(n,listt):
    for i in range(len(listt)):
        if listt[i]>=n:
            return listt[i]
    return n

def timea3DR2CDFT(n,m,p):
    bb = pyfftw.empty_aligned((n,m, p), dtype='float64')
    bf= pyfftw.empty_aligned((n,m, (p/2+1)), dtype='complex128')
    pyfftw.config.NUM_THREADS = 1 #multiprocessing.cpu_count()
    fft_object_b = pyfftw.FFTW(bb, bf,axes=(0,1,2))

    print n,m,p
    start = timer()
    fft_object_b(bb)
    end = timer()
    print end - start

#three prime numbers !      
n=3*37
m=241
p=5*19

timea3DR2CDFT(n,m,p)



# to even size :
neven=2*((n+1)/2)
meven=2*((m+1)/2)
peven=2*((p+1)/2)

timea3DR2CDFT(neven,meven,peven)


#to nearest multiple of prime
listt=listofgoodsizes()

ngood=getgoodfftwsize(n,listt)
mgood=getgoodfftwsize(m,listt)
pgood=getgoodfftwsize(p,listt)

timea3DR2CDFT(ngood,mgood,pgood)

На моем компьютере он печатает: 

111 241 95
0.180601119995
112 242 96
0.0560319423676
112 252 96
0.0564918518066
...