In [ ]:
def qs(s):
'''Compute recursive filtering scaling factor'''
return np.float64(0.00399341) + np.float64(0.4715161) * s
def dsc(d, s):
'''Rescale poles of the recursive filtering z-transform'''
q = qs(s)
r, phi = cmath.polar(d)
return cmath.rect(pow(r,np.float64(1.0)/q), phi/q)
def dsr(d, s):
'''Rescale poles in the real-axis of the recursive filtering z-transform'''
return pow(d, np.float64(1.0)/qs(s))
In [ ]:
def weights1(s):
'''Compute [van Vliet et al. 1998] 1st-order filter weights'''
s = np.float64(s)
d3 = np.float64(1.86543)
d = dsr(d3, s)
b0 = -(np.float64(1.0)-d)/d
a1 = np.float64(-1.0)/d
return b0, a1
def weights2(s):
'''Compute [van Vliet et al. 1998] 2nd-order filter weights'''
s = np.float64(s)
d1 = complex(np.float64(1.41650), np.float64(1.00829))
d = dsc(d1, s)
n2 = abs(d)
n2 *= n2
re = d.real
b0 = (np.float64(1.0)-np.float64(2.0)*re+n2)/n2
a1 = np.float64(-2.0)*re/n2
a2 = np.float64(1.0)/n2
return b0, a1, a2
def weightsk(n, k, epsilon=0.0001, theta=1.2):
'''Compute arbitrary weights'''
rho = pow(epsilon * math.sin(theta), 1./float(k*n))
x = rho * math.cos(theta)
y = rho * math.sin(theta)
b0 = 1
a1 = -2 * x
a2 = x*x + y*y
return b0, a1, a2
In [ ]:
def fwd(prol, bloc, weig):
'''Forward filter: vectors of prologue, 1D block and weights'''
for j, _ in enumerate(bloc):
bloc[j] *= weig[0]
for k, _ in enumerate(weig[1:],start=1):
if (j-k) < 0: # use data from prologue
bloc[j] -= prol[j-k]*weig[k]
else: # use data from block
bloc[j] -= bloc[j-k]*weig[k]
def FT(Prol, Bloc, weig):
'''Forward-transpose filter: matrix of prologues and block'''
FTB = np.copy(Bloc).astype(np.float64)
for row, _ in enumerate(FTB):
fwd(Prol[row], FTB[row], weig)
return FTB
def F(Prol, Bloc, weig):
'''Forward filter: matrix of prologues and block'''
return FT(Prol.transpose(), Bloc.transpose(), weig).transpose()
In [ ]:
def rev(bloc, epil, weig):
'''Reverse filter: vector of 1D block, epilogue and weights'''
for j, _ in enumerate(bloc):
j = len(bloc)-1-j
bloc[j] *= weig[0]
for k, _ in enumerate(weig[1:],start=1):
if (j+k) >= len(bloc): # use data from epilogue
bloc[j] -= epil[j+k-len(bloc)]*weig[k]
else: # use data from block
bloc[j] -= bloc[j+k]*weig[k]
def RT(Bloc, Epil, weig):
'''Reverse-transpose filter: matrix of epilogues and 2D block'''
RTB = np.copy(Bloc).astype(np.float64)
for row, _ in enumerate(RTB):
rev(RTB[row], Epil[row], weig)
return RTB
def R(Bloc, Epil, weig):
'''Reverse filter: matrix of epilogues and 2D block'''
return RT(Bloc.transpose(), Epil.transpose(), weig).transpose()
In [ ]:
def H(X, r):
'''Head operator: 2D matrix X and filter order r'''
return X[:r]
def T(X, r):
'''Tail operator: 2D matrix X and filter order r'''
return X[-r:]
def HT(X, r):
'''Head-transpose operator: 2D matrix X and filter order r'''
return X[:,:r]
def TT(X, r):
'''Tail-transpose operator: 2D matrix X and filter order r'''
return X[:,-r:]
In [ ]:
def flip(Bloc):
'''Flip both rows and columns, numpy flips give only
np.fliplr for flip_cols and np.flipud for flip_rows'''
FlipBloc = np.zeros(Bloc.shape, dtype=np.float64)
for i in range(Bloc.shape[0]):
for j in range(Bloc.shape[1]):
FlipBloc[i][j] = Bloc[Bloc.shape[0]-1-i][Bloc.shape[1]-1-j]
return FlipBloc
In [ ]:
def get_mn(X, b):
'''Get (m,n) of X given b, (m,n) are the number of blocks
inside image X of size b, X as a np.ndarray with shape as
(number of rows, number of columns)'''
m = int(math.ceil(X.shape[0] / float(b)))
n = int(math.ceil(X.shape[1] / float(b)))
return m, n
def break_blocks(X, b, m_size=None, n_size=None):
'''Break image X into blocks of size b,
where (m,n) sizes are optionally given'''
if m_size is None or n_size is None:
m_size, n_size = get_mn(X, b)
blocks = np.empty((m_size,n_size,b,b), dtype=np.float64)
for m in range(m_size):
for n in range(n_size):
blocks[m][n] = X[m*b:(m+1)*b,n*b:(n+1)*b]
return blocks
def join_blocks(blocks, b, m_size, n_size, Xshape):
'''Join blocks of size b back into returned image X
given (m,n) sizes and the shape of the original image'''
image = np.empty(Xshape, dtype=np.float64)
for m in range(m_size):
for n in range(n_size):
image[m*b:(m+1)*b,n*b:(n+1)*b] = blocks[m][n]
return image
In [ ]:
def build_alg5_matrices(b, r, wgt, width, height):
'''Build all algorithm 5 matrices and return them packed'''
Zrb = np.zeros((r, b), dtype=np.float64)
Zbr = np.zeros((b, r), dtype=np.float64)
Ir = np.identity(r, dtype=np.float64)
Ib = np.identity(b, dtype=np.float64)
Zwr = np.zeros((width, r), dtype=np.float64)
Zrh = np.zeros((r, height), dtype=np.float64)
K = np.copy(np.flipud(Ir)).astype(np.float64)
AFP = F(Ir, Zbr, wgt)
AFB = F(Zrb, Ib, wgt)
ARE = R(Zbr, Ir, wgt)
ARB = R(Ib, Zrb, wgt)
TAFP = AbF = T(AFP, r)
HARE = AbR = H(ARE, r)
HARB = H(ARB, r)
HARB_AFP = np.dot(HARB, AFP)
HARB_AFB = np.dot(HARB, AFB)
TAFPT = AbFT = TAFP.transpose()
TAFBT = T(AFB, r).transpose()
ARB_AFP = np.dot(ARB, AFP)
HARET = AbRT = HARE.transpose()
HARB_AFPT = HARB_AFP.transpose()
HARB_AFBT = HARB_AFB.transpose()
ArF = H(AFP, r)
ArR = flip(ArF)
AbarF = np.zeros((r, r))
for i in range(r):
AbarF[i][i] = wgt[0]
for j in range(i+1, r):
AbarF[j][i] = ArF[j-i-1][r-1] * wgt[0]
AbarR = flip(AbarF)
AwF = T(F(Ir, Zwr, wgt), r)
AwR = H(R(Zwr, Ir, wgt), r)
AhF = TT(FT(Ir, Zrh, wgt), r)
AhR = HT(RT(Zrh, Ir, wgt), r)
# 0, 1, 2, 3, 4, 5, 6, 7
# 8, 9, 10, 11, 12
# 13, 14, 15, 16, 17, 18, 19, 20
# 21, 22, 23, 24
return (Zrb, Zbr, Ir, Ib, AFP, AFB, AbF, AbR,\
HARB_AFP, AbFT, AbRT, HARB_AFPT, HARB_AFBT,\
ARB_AFP, TAFBT, ARE, K, ArF, ArR, AbarF, AbarR,\
AwF, AwR, AhF, AhR)
def build_alg5_carries(m_size, n_size, b, r):
'''Build all carries with zero values and return them packed'''
P = np.zeros((m_size+1,n_size,r,b), dtype=np.float64)
E = np.zeros((m_size+1,n_size,r,b), dtype=np.float64)
Pt = np.zeros((m_size,n_size+1,b,r), dtype=np.float64)
Et = np.zeros((m_size,n_size+1,b,r), dtype=np.float64)
return (P, E, Pt, Et)
In [ ]:
def alg5_stage1(m_size, n_size, r, wgt, matrices, carries, blocks):
'''Algorithm 5 stage 1'''
# unpack first two pre-computed matrices
Zrb, Zbr = matrices[:2]
P, E, Pt, Et = carries # unpack carries
for n in range(n_size):
for m in range(m_size):
B = blocks[m][n]
B = F(Zrb, B, wgt)
P[m][n] = T(B, r)
B = R(B, Zrb, wgt)
E[m][n] = H(B, r)
B = FT(Zbr, B, wgt)
Pt[m][n] = TT(B, r)
B = RT(B, Zbr, wgt)
Et[m][n] = HT(B, r)
In [ ]:
def alg5_stage23(m_size, n_size, matrices, carries):
'''Algorithm 5 stages 2 and 3'''
# unpack three pre-computed matrices
AbF, AbR, HARB_AFP = matrices[6:9]
P, E, Pt, Et = carries # unpack carries
for n in range(n_size): # In parallel for all n,
for m in range(m_size): # sequentially for each m,
P[m][n] = P[m][n] \
+ np.dot(AbF, P[m-1][n]) # eqn. (24)
for m in range(m_size-1, -1, -1): # [m_size-1,0]
E[m][n] = E[m][n] \
+ np.dot(HARB_AFP, P[m-1][n]) \
+ np.dot(AbR, E[m+1][n]) # eqn. (34)
In [ ]:
def alg5_stage45(m_size, n_size, matrices, carries):
'''Algorithm 5 stages 4 and 5'''
# unpack seven other pre-computed matrices
AbFT, AbRT, HARB_AFPT, HARB_AFBT, ARB_AFP, TAFBT, ARE = matrices[9:16]
P, E, Pt, Et = carries # unpack carries
for m in range(m_size): # In parallel for all m,
for n in range(n_size): # sequentially for each n,
Pt[m][n] = Pt[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], TAFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], TAFBT)) \
+ np.dot(Pt[m][n-1], AbFT) # eqn. (37)
for n in range(n_size-1, -1, -1): # [n_size-1,0]
Et[m][n] = Et[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], HARB_AFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], HARB_AFBT)) \
+ np.dot(Pt[m][n-1], HARB_AFPT) \
+ np.dot(Et[m][n+1], AbRT) # eqn. (39)
In [ ]:
def alg5_stage6(m_size, n_size, wgt, carries, blocks):
'''Algorithm 5 stage 6'''
P, E, Pt, Et = carries # unpack carries
for n in range(n_size):
for m in range(m_size):
B = blocks[m][n]
B = F(P[m-1][n], B, wgt)
B = R(B, E[m+1][n], wgt)
B = FT(Pt[m][n-1], B, wgt)
B = RT(B, Et[m][n+1], wgt)
blocks[m][n] = B
In [ ]:
def build_cpe_matrices(r, wgt, alg5m):
'''Build constant-padding extension matrices'''
Ir, AFP = alg5m[2], alg5m[4]
ArF, ArR, AbarF, AbarR = alg5m[17:21]
sysA, sysb = np.zeros((r*r, r*r)), np.zeros(r*r)
for i in range(r):
for j in range(r):
k = r*i + j
sysA[k][k] = 1.
for q in range(r):
for p in range(r):
u = r*p + q
sysA[k][u] -= ArR[j][q] * ArF[p][i]
sysb[k] = AbarR[j][i]
sysx = linalg.lu_solve(linalg.lu_factor(sysA), sysb)
SRF = np.zeros((r, r))
for i in range(r):
for j in range(r):
k = r*i + j
SRF[j][i] = sysx[k]
IArF, IArR = Ir - ArF, Ir - ArR
SF, SR = np.linalg.inv(IArF), np.linalg.inv(IArR)
SFAbarF, SRFArF = np.dot(SF, AbarF), np.dot(SRF, ArF)
SRAbaRSRFArFSFAbarF = np.dot(np.dot((np.dot(SR, AbarR) - np.dot(SRF, ArF)), SF), AbarF)
AbarFSFT, ArFSRFT = SFAbarF.transpose(), SRFArF.transpose()
AbarFSFSRAbaRSRFArFT = SRAbaRSRFArFSFAbarF.transpose()
return (SFAbarF, SRFArF, SRAbaRSRFArFSFAbarF, AbarFSFT, ArFSRFT, AbarFSFSRAbaRSRFArFT)
In [ ]:
def alg5_cpe_stage1(m_size, n_size, r, wgt, matrices, carries, blocks):
'''Algorithm 5 constant-padding extension stage 1
The algorithm stage 5.1 for the CPE differs from the original 5.1
for zero-padding by saving in outside carries the X and Zhat values
needed by the next fixing stages'''
Zrb, Zbr = matrices[:2] # unpack first two pre-computed matrices
P, E, Pt, Et = carries # unpack carries
for n in range(n_size):
for m in range(m_size):
B = blocks[m][n]
if m == 0:
P[-1][n] = np.tile(H(B, 1), (r, 1)) # replicate first row r times
elif m == m_size-1:
E[m_size][n] = np.tile(T(B, 1), (r, 1)) # replicate last row r times
B = F(Zrb, B, wgt)
P[m][n] = T(B, r)
B = R(B, Zrb, wgt)
E[m][n] = H(B, r)
if n == 0:
Pt[m][-1] = np.tile(HT(B, 1), (1, r)) # replicate first column r times
elif n == n_size-1:
Et[m][n_size] = np.tile(TT(B, 1), (1, r)) # replicate last column r times
B = FT(Zbr, B, wgt)
Pt[m][n] = TT(B, r)
B = RT(B, Zbr, wgt)
Et[m][n] = HT(B, r)
In [ ]:
def alg5_cpe_stage23(m_size, n_size, matrices, cpe_matrices, carries):
'''Algorithm 5 constant-padding extension stages 2 and 3'''
AbF, AbR, HARB_AFP = matrices[6:9] # unpack three pre-computed matrices
SFAbarF, SRFArF, SRAbaRSRFArFSFAbarF = cpe_matrices[:3] # unpack three pre-computed CPE matrices
P, E, Pt, Et = carries # unpack carries
for n in range(n_size): # In parallel for all n,
# prior to enter the image from top to bottom (a.k.a. m = -1 index)
P[-1][n] = np.dot(SFAbarF, P[-1][n]) # *NEW PAPER* eqn. (13)
for m in range(m_size): # sequentially for each m,
P[m][n] = P[m][n] \
+ np.dot(AbF, P[m-1][n]) # eqn. (24)
# prior to enter the image from bottom to top (a.k.a. m = m_size index)
E[m_size][n] = np.dot(SRFArF, P[m_size-1][n]) \
+ np.dot(SRAbaRSRFArFSFAbarF, E[m_size][n]) # *NEW PAPER* eqn. (22)
for m in range(m_size-1, -1, -1): # [m_size-1,0]
E[m][n] = E[m][n] \
+ np.dot(HARB_AFP, P[m-1][n]) \
+ np.dot(AbR, E[m+1][n]) # eqn. (34)
In [ ]:
def alg5_cpe_stage45(m_size, n_size, r, matrices, cpe_matrices, carries):
'''Algorithm 5 constant-padding extension stages 4 and 5'''
# unpack seven other pre-computed matrices
AbFT, AbRT, HARB_AFPT, HARB_AFBT, ARB_AFP, TAFBT, ARE = matrices[9:16]
# unpack three pre-computed CPE matrices
AbarFSFT, ArFSRFT, AbarFSFSRAbaRSRFArFT = cpe_matrices[3:]
P, E, Pt, Et = carries # unpack carries
for m in range(m_size): # In parallel for all m,
# prior to enter the image from left to right (a.k.a. n = -1 index)
# compute one-up corner north west (cnw) from fixed prologues at the first column of blocks (n = 0)
cnw = np.tile(HT(P[m-1][0], 1), (1, r)) # replicate first column r times
# compute one-down corner south west (csw) from fixed epilogues at the first column of blocks (n = 0)
csw = np.tile(HT(E[m+1][0], 1), (1, r)) # replicate first column r times
# fixing Pt_-1(Zhat) -> Pt_-1(Z)
Pt[m][-1] = Pt[m][-1] \
+ np.dot(ARB_AFP, cnw) \
+ np.dot(ARE, csw) # *SIMPLIFIED* *OLD* eqn. (37)
# fixing Pt_-1(Z) -> Pt_-1(U)
Pt[m][-1] = np.dot(Pt[m][-1], AbarFSFT) # *NEW PAPER* *TRANSPOSED* eqn. (13)
for n in range(n_size): # sequentially for each n,
Pt[m][n] = Pt[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], TAFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], TAFBT)) \
+ np.dot(Pt[m][n-1], AbFT) # eqn. (37)
# prior to enter the image from right to left (a.k.a. n = n_size index)
# compute one-up corner north east (cne) from fixed prologues at the last column of blocks (n = n_size-1)
cne = np.tile(TT(P[m-1][n_size-1], 1), (1, r)) # replicate last column r times
# compute one-down corner south east (cse) from fixed epilogues at the last column of blocks (n = n_size-1)
cse = np.tile(TT(E[m+1][n_size-1], 1), (1, r)) # replicate last column r times
# fixing Et_N(Zhat) -> Et_N(Z)
Et[m][n_size] = Et[m][n_size] \
+ np.dot(ARB_AFP, cne) \
+ np.dot(ARE, cse) # *SIMPLIFIED* *OLD* eqn. (37)
# fixing Et_N(Z) -> Et_N(U)
Et[m][n_size] = np.dot(Pt[m][n_size-1], ArFSRFT) \
+ np.dot(Et[m][n_size], AbarFSFSRAbaRSRFArFT) # *NEW PAPER* *TRANSPOSED* eqn. (22)
for n in range(n_size-1, -1, -1): # [n_size-1,0]
Et[m][n] = Et[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], HARB_AFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], HARB_AFBT)) \
+ np.dot(Pt[m][n-1], HARB_AFPT) \
+ np.dot(Et[m][n+1], AbRT) # eqn. (39)
In [ ]:
def build_pe_matrices(r, wgt, alg5m):
'''Build periodic extension matrices'''
Ir = alg5m[2]
AwF, AwR, AhF, AhR = alg5m[21:25]
IAwF = np.linalg.inv(Ir - AwF)
IAwR = np.linalg.inv(Ir - AwR)
IAhF = np.linalg.inv(Ir - AhF).transpose()
IAhR = np.linalg.inv(Ir - AhR).transpose()
return (IAhF, IAhR, IAwF, IAwR)
In [ ]:
def alg5_pe_stage23(m_size, n_size, matrices, pe_matrices, carries):
'''Algorithm 5 periodic extension stages 2 and 3'''
Zrb = matrices[0]
AbF, AbR, HARB_AFP = matrices[6:9] # unpack three pre-computed matrices
IAhF, IAhR = pe_matrices[:2] # unpack first two pre-computed PE matrices
P, E, Pt, Et = carries # unpack carries
for n in range(n_size): # In parallel for all n,
PM1Y = np.copy(Zrb)
for m in range(m_size): # sequentially for each m,
PM1Y = P[m][n] \
+ np.dot(AbF, PM1Y) # eqn. (24)
P[-1][n] = np.dot(IAhF, PM1Y) # *NEW PAPER* eqn. (39)
for m in range(m_size): # sequentially for each m,
P[m][n] = P[m][n] \
+ np.dot(AbF, P[m-1][n]) # eqn. (24)
E0Z = np.copy(Zrb)
for m in range(m_size-1, -1, -1): # [m_size-1,0]
E0Z = E[m][n] \
+ np.dot(HARB_AFP, P[m-1][n]) \
+ np.dot(AbR, E0Z) # eqn. (34)
E[m_size][n] = np.dot(IAhR, E0Z) # *NEW PAPER* eqn. (41)
for m in range(m_size-1, -1, -1): # [m_size-1,0]
E[m][n] = E[m][n] \
+ np.dot(HARB_AFP, P[m-1][n]) \
+ np.dot(AbR, E[m+1][n]) # eqn. (34)
In [ ]:
def alg5_pe_stage45(m_size, n_size, r, matrices, pe_matrices, carries):
'''Algorithm 5 periodic extension stages 4 and 5'''
# unpack seven other pre-computed matrices
AbFT, AbRT, HARB_AFPT, HARB_AFBT, ARB_AFP, TAFBT, ARE = matrices[9:16]
Zbr = matrices[1]
IAwF, IAwR = pe_matrices[2:] # unpack second two pre-computed PE matrices
P, E, Pt, Et = carries # unpack carries
for m in range(m_size): # In parallel for all m,
PtNU = np.copy(Zbr)
for n in range(n_size): # sequentially for each n,
PtNU = Pt[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], TAFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], TAFBT)) \
+ np.dot(PtNU, AbFT) # eqn. (37)
Pt[m][-1] = np.dot(PtNU, IAwF) # *NEW PAPER* eqn. (39)
for n in range(n_size): # sequentially for each n,
Pt[m][n] = Pt[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], TAFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], TAFBT)) \
+ np.dot(Pt[m][n-1], AbFT) # eqn. (37)
Et0V = np.copy(Zbr)
for n in range(n_size-1, -1, -1): # [n_size-1,0]
Et0V = Et[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], HARB_AFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], HARB_AFBT)) \
+ np.dot(Pt[m][n-1], HARB_AFPT) \
+ np.dot(Et0V, AbRT) # eqn. (39)
Et[m][n_size] = np.dot(Et0V, IAwR) # *NEW PAPER* eqn. (41)
for n in range(n_size-1, -1, -1): # [n_size-1,0]
Et[m][n] = Et[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], HARB_AFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], HARB_AFBT)) \
+ np.dot(Pt[m][n-1], HARB_AFPT) \
+ np.dot(Et[m][n+1], AbRT) # eqn. (39)
In [ ]:
def build_epe_matrices(r, wgt, alg5m):
'''Build even-periodic extension matrices'''
Ir, AbF, AbR = alg5m[2], alg5m[6], alg5m[7]
K, ArF, ArR, AbarF, AbarR, AwF, AwR, AhF, AhR = alg5m[16:25]
L = np.dot(np.linalg.inv(K - ArR), AbarR)
IA2wF = np.linalg.inv(Ir - np.dot(AwF, AwF))
IA2hF = np.linalg.inv(Ir - np.dot(AhF, AhF))
AbarFi = np.linalg.inv(AbarF)
M1w = np.dot(IA2wF, np.dot(K, np.dot(AbarFi, (Ir - np.dot(ArF, ArR)))))
M2w = np.dot(IA2wF, (AwF + np.dot(K, np.dot(AwR, np.dot(AbarFi, np.dot(ArF, AbarR))))))
M1h = np.dot(IA2hF, np.dot(K, np.dot(AbarFi, (Ir - np.dot(ArF, ArR)))))
M2h = np.dot(IA2hF, (AhF + np.dot(K, np.dot(AhR, np.dot(AbarFi, np.dot(ArF, AbarR))))))
return L, M1w, M2w, M1h, M2h
In [ ]:
def alg5_epe_stage23(m_size, n_size, matrices, epe_matrices, carries):
'''Algorithm 5 even-periodic extension stages 2 and 3'''
# unpack five pre-computed matrices
Zrb, AhF = matrices[0], matrices[23]
AbF, AbR, HARB_AFP = matrices[6:9]
# unpack three pre-computed EPE matrices
L, M1h, M2h = epe_matrices[0], epe_matrices[3], epe_matrices[4]
P, E, Pt, Et = carries # unpack carries
tP = np.copy(P) # temp P
for n in range(n_size): # In parallel for all n,
E0Z = np.copy(Zrb)
for m in range(m_size): # sequentially for each m,
tP[m][n] = tP[m][n] \
+ np.dot(AbF, tP[m-1][n]) # eqn. (24)
for m in range(m_size-1, -1, -1): # [m_size-1,0]
E0Z = E[m][n] \
+ np.dot(HARB_AFP, tP[m-1][n]) \
+ np.dot(AbR, E0Z) # eqn. (34)
P[-1][n] = np.dot(M1h, E0Z) \
+ np.dot(M2h, tP[m_size-1][n]) # *NEW PAPER* eqn. (57)
E[m_size][n] = np.dot(L, \
(np.dot(AhF, P[-1][n]) + tP[m_size-1][n])) # *NEW PAPER* eqn. (62)
for m in range(m_size): # sequentially for each m,
P[m][n] = P[m][n] \
+ np.dot(AbF, P[m-1][n]) # eqn. (24)
for m in range(m_size-1, -1, -1): # [m_size-1,0]
E[m][n] = E[m][n] \
+ np.dot(HARB_AFP, P[m-1][n]) \
+ np.dot(AbR, E[m+1][n]) # eqn. (34)
In [ ]:
def alg5_epe_stage45(m_size, n_size, r, matrices, epe_matrices, carries):
'''Algorithm 5 even-periodic extension stages 4 and 5'''
# unpack nine other pre-computed matrices
Zbr, AwF = matrices[1], matrices[21]
AbFT, AbRT, HARB_AFPT, HARB_AFBT, ARB_AFP, TAFBT, ARE = matrices[9:16]
# unpack three pre-computed EPE matrices
L, M1w, M2w = epe_matrices[0:3]
P, E, Pt, Et = carries # unpack carries
tPt = np.copy(Pt) # temp Pt
for m in range(m_size): # In parallel for all m,
Et0V = np.copy(Zbr)
for n in range(n_size): # sequentially for each n,
tPt[m][n] = tPt[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], TAFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], TAFBT)) \
+ np.dot(tPt[m][n-1], AbFT) # eqn. (37)
for n in range(n_size-1, -1, -1): # [n_size-1,0]
Et0V = Et[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], HARB_AFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], HARB_AFBT)) \
+ np.dot(tPt[m][n-1], HARB_AFPT) \
+ np.dot(Et0V, AbRT) # eqn. (39)
Pt[m][-1] = np.dot(Et0V, M1w.transpose()) \
+ np.dot(tPt[m][n_size-1], M2w.transpose()) # *NEW PAPER* eqn. (57)
Et[m][n_size] = np.dot((np.dot(Pt[m][-1], \
AwF.transpose()) + tPt[m][n_size-1]), L.transpose()) # *NEW PAPER* eqn. (62)
for n in range(n_size): # sequentially for each n,
Pt[m][n] = Pt[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], TAFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], TAFBT)) \
+ np.dot(Pt[m][n-1], AbFT) # eqn. (37)
for n in range(n_size-1, -1, -1): # [n_size-1,0]
Et[m][n] = Et[m][n] \
+ np.dot(ARB_AFP, np.dot(P[m-1][n], HARB_AFBT)) \
+ np.dot(ARE, np.dot(E[m+1][n], HARB_AFBT)) \
+ np.dot(Pt[m][n-1], HARB_AFPT) \
+ np.dot(Et[m][n+1], AbRT) # eqn. (39)
In [ ]:
print 'All functions loaded successfully!'