Last active
August 13, 2018 13:34
-
-
Save Foadsf/80abba07ea892e50da3df7e318844f33 to your computer and use it in GitHub Desktop.
Cauchy product of multivariate finite power series (polynomials) represented as discrete convolution of NumPy ndarrays
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
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 xpndzr(A,D): | |
D1=D-A.shape | |
T=np.zeros((D.shape[0],2),dtype=D.dtype) | |
T[:,1]=D1[:] | |
return np.pad(A, T, 'constant', constant_values=0) | |
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))))) | |
D1=np.array([4,5]) | |
D2=np.array([2,3]) | |
A=np.random.randint(10,size=D1) | |
B=np.random.randint(10,size=D2) | |
K=np.array([4,6])+1 | |
print(conv(A,B,K)) |
I am trying to parse your formulas and your code, and am trying to come up with an example implementation as simple as possible. I could come up with:
import numpy as np
def crop(A, start=0, stop=-1):
"""
crop n-dimensional array along all dimensions:
A -> A[start:stop, start:stop, ..., start:stop]
"""
return A[[slice(start, stop)] * A.ndim]
def flipall(A):
"""
flip n-dimensional array along all dimensions:
A -> A[::-1, ::-1, ..., ::-1]
"""
return A[[slice(None, None, -1)] * A.ndim]
def conv(A, B, j):
return np.sum(np.conj(crop(A, stop=j)) * flipall(crop(B, stop=j)))
A=np.random.randint(10,size=(4, 4))
B=np.random.randint(10,size=(4, 4))
print(A)
print(B)
print(conv(A, B, 1))
print(conv(A, B, 2))
does this do what you want?
@nils-werner I edited the original post and changed J
to K
to be consistent with the original notation. K=[k1,...,kn]
where for all 0<j<=n
we have 0<=kj<=oj
@nils-werner Thanks alot for the sample code, it is awesome. just some questions:
- the
crop
function which I have explained here should accept two 1D arraysD1=[d11,...,d1n]
andD2=[d21,...,d2n]
of non-negative integers, for the last two arguments. whereA.ndim=D1.shape[0]=D2.shape[0]
and for all0<j<=n
we should have0<=d1j<=d2j<=A.shape[j]
. What I do not understand is how the default values in your code are just integers. shouldn't they be ndarrays? - I'm not sure your
conv
function is same as mine. I'm going to test it now. Formerj
orJ
and nowK
, is an array with the same shape asD1
andD2
. but in your code it is an integer. A
andB
are not necessarily square/cubical. they can have any arbitrary shapesD1
andD2
.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
What's
J
?