Как реализовать многомерную версию NymPy?
После этого сообщения мне посоветовали задать другой вопрос на основе MCVE. Моя цель - реализовать свертку NumPy для произвольно сформированных входных массивов. Учтите, что я использую термин свертка в контексте произведения Коши многомерных степенных рядов (многопараметрическое многочленное умножение). Функции SciPy, такие как signal.convolve, ndimage.convolve или ndimage.filters.convolve, не работают для меня, как я объяснил здесь.
Рассмотрим два неквадратных 2D массива NumPy A и B:
D1=np.array([4,5])
D2=np.array([2,3])
A=np.random.randint(10,size=D1)
B=np.random.randint(10,size=D2)
, например:
[[1 4 4 2 7]
[6 1 7 5 3]
[1 4 3 4 8]
[7 5 8 3 3]]
[[2 2 3]
[5 2 9]]
Теперь я способный вычислять элементы C=convolve(A,B) с помощью conv(A,B,K):
def crop(A,D1,D2):
return A[tuple(slice(D1[i], D2[i]) for i in range(A.ndim))]
def sumall(A):
sum1=A
for k in range(A.ndim):
sum1 = np.sum(sum1,axis=0)
return sum1
def flipall(A):
return A[[slice(None, None, -1)] * A.ndim]
def conv(A,B,K):
D0=np.zeros(K.shape,dtype=K.dtype)
return sumall(np.multiply(crop(A,np.maximum(D0,np.minimum(A.shape,K-B.shape)) \
,np.minimum(A.shape,K)), \
flipall(crop(B,np.maximum(D0,np.minimum(B.shape,K-A.shape)) \
,np.minimum(B.shape,K)))))
Пример Fow для K=np.array([0,0])+1, conve(A,B,K) результатов 1*2=2, для K=np.array([1,0])+1 результатов 5*1+2*6=17, для K=np.array([0,1])+1 есть 2*4+1*2=10, а для K=np.array([1,1])+1 дает 4*5+6*2+1*1+1*2=36:
[[2 10 ...]
[17 36 ...]
... ]]
теперь, если бы я знал размерность A и B, я мог бы вложить некоторые из циклов для заполнения C, но это не относится ко мне. Как я могу использовать функцию conv для заполнения C ndarray с формой C.shape=A.shape+B.shape-1 без использования для циклов?
1 ответ
Сначала нам нужны данные для сравнения. Итак, давайте разберем пространство ввода:
d0, d1 = np.array(A.shape) + np.array(B.shape) - 1
input_space = np.array(np.meshgrid(np.arange(d0), np.arange(d1))).T.reshape(-1, 2)
# array([[0, 0],
# [0, 1],
# [0, 2],
# [0, 3],
# [0, 4],
# [0, 5],
# [0, 6],
# [1, 0],
# [1, 1],
# ...
# [4, 5],
# [4, 6]])
и вычислим вашу свертку над этим пространством:
out = np.zeros((d0, d1))
for K in input_space:
out[tuple(K)] = conv(A, B, K + 1)
out
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
Хорошо, теперь, когда мы знаем, какие значения ожидать, давайте посмотрим, мы можем получить scipy и astropy, чтобы дать нам те же значения:
import scipy.signal
scipy.signal.convolve2d(A, B) # only 2D!
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
import astropy.convolution
astropy.convolution.convolve_fft(
np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
B,
normalize_kernel=False
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
astropy.convolution.convolve(
np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
np.pad(B, pad_width=((0, 1), (0, 0)), mode='constant'),
normalize_kernel=False
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
import scipy
scipy.ndimage.filters.convolve(
np.pad(A, pad_width=((0, 1), (0, 2)), mode='constant'),
B,
mode='constant',
cval=0.0,
origin=-1
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
scipy.ndimage.filters.convolve(
np.pad(A, pad_width=((1, 0), (1, 1)), mode='constant'),
B,
mode='constant',
cval=0.0
)
# array([[ 2., 10., 19., 24., 30., 20., 21.],
# [ 17., 36., 71., 81., 112., 53., 72.],
# [ 32., 27., 108., 74., 121., 79., 51.],
# [ 19., 46., 79., 99., 111., 67., 81.],
# [ 35., 39., 113., 76., 93., 33., 27.]])
Как вы можете видеть, это просто вопрос правильной нормализации и дополнения, и вы можете просто использовать любая из этих библиотек.
Я рекомендую использовать astropy.convolution.convolve_fft, поскольку он (будучи основанным на FFT), вероятно, является самым быстрым.
-
1это круто. Вы ответили на мой первоначальный вопрос. Я собираюсь попробовать их и прийти сюда. Я ценю вашу поддержку. – Foad 13 August 2018 в 16:03
-
2Теперь я понимаю, что вы сделали. функции, приведенные выше, фактически предназначены для обработки изображений, поэтому вы расширили их с помощью numpy.pad, а затем использовали функции scipy или astropy. но для этого требуется знать A и B. если мы автоматизируем эту часть, тогда я получаю то, что мне нужно. – Foad 13 August 2018 в 16:25
-
3Моя конечная цель - умножить многомерные конечные степенные ряды. в какой-то момент мне нужно увеличить степень сери для конкретной переменной. Хотя делать это вручную возможно, это будет неэффективно. Было бы идеально построить функцию Convolve, чтобы взять A и B и дать результат. – Foad 13 August 2018 в 16:33
-
4Вот что делают все эти библиотеки! Вы даете им A и B (возможно, набиваете некоторые нули здесь и там), и вы получаете результат. – Nils Werner 13 August 2018 в 16:33