shuren/fdfd_tools
1
0
Fork 0
forked from jan/fdfd_tools
Code Issues Pull Requests Releases Wiki Activity

cleanup

Browse Source
This commit is contained in:
jan 2018-01-08 23:28:57 -08:00
parent 4067766478
commit c4cbdff751
1 changed files with 63 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

125
fdfd_tools/bloch.py
View File

@ -447,15 +447,10 @@ def eigsolve(num_modes: int,
continue continue
break break
def rtrace_AtB(A, B):
return real(numpy.sum(A.conj() * B))
def symmetrize(A):
return (A + A.conj().T) * 0.5
max_iters = 10000 max_iters = 10000
for iter in range(max_iters): for iter in range(max_iters):
U = numpy.linalg.inv(Z.conj().T @ Z) ZtZ = Z.conj().T @ Z
U = numpy.linalg.inv(ZtZ)
AZ = scipy_op @ Z AZ = scipy_op @ Z
AZU = AZ @ U AZU = AZ @ U
ZtAZU = Z.conj().T @ AZU ZtAZU = Z.conj().T @ AZU
@ -469,47 +464,44 @@ def eigsolve(num_modes: int,
break break
KG = scipy_iop @ G KG = scipy_iop @ G
traceGtKG = rtrace_AtB(G, KG) traceGtKG = _rtrace_AtB(G, KG)
gamma_numerator = traceGtKG
reset_iters = 100 reset_iters = 100 # TODO
if prev_traceGtKG == 0 or iter % reset_iters == 0: if prev_traceGtKG == 0 or iter % reset_iters == 0:
print('RESET!') logger.inf('CG reset')
gamma = 0 gamma = 0
else: else:
gamma = gamma_numerator / prev_traceGtKG gamma = traceGtKG / prev_traceGtKG
D = gamma * d_scale * D + KG D = gamma * d_scale * D + KG
d_scale = numpy.sqrt(rtrace_AtB(D, D)) / num_modes d_scale = numpy.sqrt(_rtrace_AtB(D, D)) / num_modes
D /= d_scale D /= d_scale
ZtAZ = Z.conj().T @ AZ
AD = scipy_op @ D AD = scipy_op @ D
DtD = D.conj().T @ D DtD = D.conj().T @ D
DtAD = D.conj().T @ AD DtAD = D.conj().T @ AD
ZtD = Z.conj().T @ D symZtD = _symmetrize(Z.conj().T @ D)
ZtAD = Z.conj().T @ AD symZtAD = _symmetrize(Z.conj().T @ AD)
symZtD = symmetrize(ZtD)
symZtAD = symmetrize(ZtAD)
'''
U_sZtD = U @ symZtD U_sZtD = U @ symZtD
dE = 2.0 * (rtrace_AtB(U, symZtAD) - rtrace_AtB(ZtAZU, U_sZtD)) dE = 2.0 * (_rtrace_AtB(U, symZtAD) -
_rtrace_AtB(ZtAZU, U_sZtD))
S2 = DtD - 4 * symZtD @ U_sZtD d2E = 2 * (_rtrace_AtB(U, DtAD) -
d2E = 2 * (rtrace_AtB(U, DtAD) - _rtrace_AtB(ZtAZU, U @ (DtD - 4 * symZtD @ U_sZtD)) -
rtrace_AtB(ZtAZU, U @ S2) - 4 * _rtrace_AtB(U, symZtAD @ U_sZtD))
4 * rtrace_AtB(U, symZtAD @ U_sZtD))
# Newton-Raphson to find a root of the first derivative: # Newton-Raphson to find a root of the first derivative:
theta = -dE/d2E theta = -dE/d2E
if d2E < 0 or abs(theta) >= pi: if d2E < 0 or abs(theta) >= pi:
theta = -abs(prev_theta) * numpy.sign(dE) theta = -abs(prev_theta) * numpy.sign(dE)
'''
# ZtAZU * ZtZ = ZtAZ for use in line search
ZtZ = Z.conj().T @ Z
ZtAZ = ZtAZU @ ZtZ.conj().T
def Qi_func(theta, memo=[None, None]): def Qi_func(theta, memo=[None, None]):
if memo[0] == theta: if memo[0] == theta:
@ -525,10 +517,10 @@ def eigsolve(num_modes: int,
# if c or s small, taylor expand # if c or s small, taylor expand
if c < 1e-4 * s and c != 0: if c < 1e-4 * s and c != 0:
Qi = numpy.linalg.inv(DtD) Qi = numpy.linalg.inv(DtD)
Qi = Qi / (s*s) - 2*c/(s*s*s) * (Qi @ symZtD.conj().T @ Qi.conj().T) Qi = Qi / (s*s) - 2*c/(s*s*s) * (Qi @ (Qi @ symZtD).conj().T)
elif s < 1e-4 * c and s != 0: elif s < 1e-4 * c and s != 0:
Qi = numpy.linalg.inv(ZtZ) Qi = numpy.linalg.inv(ZtZ)
Qi = Qi / (c*c) - 2*s/(c*c*c) * (Qi @ symZtD.conj().T @ Qi.conj().T) Qi = Qi / (c*c) - 2*s/(c*c*c) * (Qi @ (Qi @ symZtD).conj().T)
else: else:
raise Exception('Inexplicable singularity in trace_func') raise Exception('Inexplicable singularity in trace_func')
memo[0] = theta memo[0] = theta
@ -540,22 +532,24 @@ def eigsolve(num_modes: int,
s = numpy.sin(theta) s = numpy.sin(theta)
Qi = Qi_func(theta) Qi = Qi_func(theta)
R = c*c * ZtAZ + s*s * DtAD + 2*s*c * symZtAD R = c*c * ZtAZ + s*s * DtAD + 2*s*c * symZtAD
trace = rtrace_AtB(R, Qi) trace = _rtrace_AtB(R, Qi)
return numpy.abs(trace) return numpy.abs(trace)
#def trace_deriv(theta): '''
# Qi = Qi_func(theta) def trace_deriv(theta):
# c2 = numpy.cos(2 * theta) Qi = Qi_func(theta)
# s2 = numpy.sin(2 * theta) c2 = numpy.cos(2 * theta)
# F = -0.5*s2 * (ZtAZ - DtAD) + c2 * symZtAD s2 = numpy.sin(2 * theta)
# trace_deriv = rtrace_AtB(Qi, F) F = -0.5*s2 * (ZtAZ - DtAD) + c2 * symZtAD
trace_deriv = _rtrace_AtB(Qi, F)
# G = Qi @ F.conj().T @ Qi.conj().T G = Qi @ F.conj().T @ Qi.conj().T
# H = -0.5*s2 * (ZtZ - DtD) + c2 * symZtD H = -0.5*s2 * (ZtZ - DtD) + c2 * symZtD
# trace_deriv -= rtrace_AtB(G, H) trace_deriv -= _rtrace_AtB(G, H)
# trace_deriv *= 2 trace_deriv *= 2
# return trace_deriv * sgn return trace_deriv * sgn
'''
''' '''
theta, new_E, new_dE = linmin(theta, E, dE, 0.1, min(tolerance, 1e-6), 1e-14, 0, -numpy.sign(dE) * K_PI, trace_func) theta, new_E, new_dE = linmin(theta, E, dE, 0.1, min(tolerance, 1e-6), 1e-14, 0, -numpy.sign(dE) * K_PI, trace_func)
@ -597,29 +591,36 @@ def eigsolve(num_modes: int,
order = numpy.argsort(numpy.abs(eigvals)) order = numpy.argsort(numpy.abs(eigvals))
return eigvals[order], eigvecs.T[order] return eigvals[order], eigvecs.T[order]
#def linmin(x_guess, f0, df0, x_max, f_tol=0.1, df_tol=min(tolerance, 1e-6), x_tol=1e-14, x_min=0, linmin_func): #def linmin(x_guess, f0, df0, x_max, f_tol=0.1, df_tol=min(tolerance, 1e-6), x_tol=1e-14, x_min=0, linmin_func):
# if df0 > 0: # if df0 > 0:
# x0, f0, df0 = linmin(-x_guess, f0, -df0, -x_max, f_tol, df_tol, x_tol, -x_min, lambda q, dq: -linmin_func(q, dq)) # x0, f0, df0 = linmin(-x_guess, f0, -df0, -x_max, f_tol, df_tol, x_tol, -x_min, lambda q, dq: -linmin_func(q, dq))
# return -x0, f0, -df0 # return -x0, f0, -df0
# elif df0 == 0: # elif df0 == 0:
# return 0, f0, df0 # return 0, f0, df0
# else: # else:
# x = x_guess # x = x_guess
# fx = f0 # fx = f0
# dfx = df0 # dfx = df0
# isave = numpy.zeros((2,), numpy.intc) # isave = numpy.zeros((2,), numpy.intc)
# dsave = numpy.zeros((13,), float) # dsave = numpy.zeros((13,), float)
# x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task, # x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
# x_min, x_max, isave, dsave) # x_min, x_max, isave, dsave)
# for i in range(int(1e6)): # for i in range(int(1e6)):
# if task != 'F': # if task != 'F':
# logging.info('search converged in {} iterations'.format(i)) # logging.info('search converged in {} iterations'.format(i))
# break # break
# fx = f(x, dfx) # fx = f(x, dfx)
# x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task, # x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
# x_min, x_max, isave, dsave) # x_min, x_max, isave, dsave)
# return x, fx, dfx # return x, fx, dfx
def _rtrace_AtB(A, B):
return real(numpy.sum(A.conj() * B))
def _symmetrize(A):
return (A + A.conj().T) * 0.5