mirror of
https://github.com/vale981/stocproc
synced 2025-03-05 09:41:42 -05:00
work in progress, figuring out the optimal interpolation strategy
This commit is contained in:
Richard Hartmann
parent
1fcb195aa7
commit
ddd99b4272
4 changed files with 677 additions and 78 deletions
|
|
@ -1,11 +1,16 @@
|
|||
import numpy as np
|
||||
from scipy.linalg import eigh as scipy_eigh
|
||||
from scipy.special import gamma
|
||||
from scipy.optimize import minimize
|
||||
import time
|
||||
from . import stocproc_c
|
||||
from . import gquad
|
||||
from . import tools
|
||||
|
||||
import logging
|
||||
log = logging.getLogger(__name__)
|
||||
|
||||
def solve_hom_fredholm(r, w, eig_val_min):
|
||||
def solve_hom_fredholm(r, w, eig_val_min=None):
|
||||
r"""Solves the discrete homogeneous Fredholm equation of the second kind
|
||||
|
||||
.. math:: \int_0^{t_\mathrm{max}} \mathrm{d}s R(t-s) u(s) = \lambda u(t)
|
||||
|
|
@ -41,17 +46,52 @@ def solve_hom_fredholm(r, w, eig_val_min):
|
|||
w_sqrt = np.sqrt(w)
|
||||
r = w_sqrt.reshape(n,1) * r * w_sqrt.reshape(1,n)
|
||||
eig_val, eig_vec = scipy_eigh(r, overwrite_a=True) # eig_vals in ascending
|
||||
min_idx = sum(eig_val < eig_val_min)
|
||||
|
||||
if eig_val_min is None:
|
||||
min_idx = 0
|
||||
else:
|
||||
min_idx = sum(eig_val < eig_val_min)
|
||||
|
||||
eig_val = eig_val[min_idx:][::-1]
|
||||
eig_vec = eig_vec[:, min_idx:][:, ::-1]
|
||||
log.debug("discrete fredholm equation of size {} solved [{:.2e}]".format(n, time.time()-t0))
|
||||
|
||||
num_of_functions = len(eig_val)
|
||||
log.debug("use {} / {} eigenfunctions (sig_min = {})".format(num_of_functions, len(w), np.sqrt(eig_val_min)))
|
||||
log.debug("use {} / {} eigenfunctions".format(len(eig_val), len(w)))
|
||||
eig_vec = np.reshape(1/w_sqrt, (n,1)) * eig_vec
|
||||
|
||||
log.debug("discrete fredholm equation of size {} solved [{:.2e}]".format(n, time.time() - t0))
|
||||
return eig_val, eig_vec
|
||||
|
||||
def align_eig_vec(eig_vec):
|
||||
for i in range(eig_vec.shape[1]):
|
||||
phase = np.exp(1j * np.arctan2(np.real(eig_vec[0,i]), np.imag(eig_vec[0,i])))
|
||||
eig_vec[:, i] /= phase
|
||||
|
||||
def opt_fredholm(kernel, lam, t, u, meth, ntimes):
|
||||
def fopt(x, p, k):
|
||||
n = len(x)//2 + 1
|
||||
lam = x[0]
|
||||
u = x[1:n] + 1j*x[n:]
|
||||
return np.log10((np.sum(np.abs(k.dot(u) - lam * u) ** p)) ** (1 / p))
|
||||
|
||||
for i in range(ntimes):
|
||||
print("iter {}/{} lam {}".format(i+1, ntimes, lam))
|
||||
u_intp = tools.ComplexInterpolatedUnivariateSpline(t, u)
|
||||
ng = 2*(len(t)-1)+1
|
||||
t, w = meth(t[-1], ng)
|
||||
u_new = u_intp(t)
|
||||
k_new = kernel(t.reshape(-1,1)-t.reshape(1,-1))*w.reshape(-1,1)
|
||||
p = 5
|
||||
r = minimize(fopt, x0=np.hstack((lam, u_new.real, u_new.imag)), args=(p, k_new), method='CG')
|
||||
x = r.x
|
||||
n = len(x) // 2 + 1
|
||||
lam = x[0]
|
||||
u = x[1:n] + 1j * x[n:]
|
||||
print("func val {} ng {}".format(r.fun, ng))
|
||||
|
||||
return lam, u, t
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def get_mid_point_weights_times(t_max, num_grid_points):
|
||||
r"""Returns the time gridpoints and wiehgts for numeric integration via **mid point rule**.
|
||||
|
||||
|
|
@ -85,8 +125,8 @@ def get_trapezoidal_weights_times(t_max, num_grid_points):
|
|||
return t, w
|
||||
|
||||
def get_simpson_weights_times(t_max, num_grid_points):
|
||||
if num_grid_points % 2 == 0:
|
||||
raise RuntimeError("simpson weight need odd number of grid points, but git ng={}".format(num_grid_points))
|
||||
if num_grid_points % 2 != 1:
|
||||
raise RuntimeError("simpson weights needs grid points ng such that ng = 2*k+1, but git ng={}".format(num_grid_points))
|
||||
# generate mid points
|
||||
t, delta_t = np.linspace(0, t_max, num_grid_points, retstep = True)
|
||||
# equal weights for grid points
|
||||
|
|
@ -95,4 +135,93 @@ def get_simpson_weights_times(t_max, num_grid_points):
|
|||
w[1::2] = 4/3*delta_t
|
||||
w[0] = 1/3*delta_t
|
||||
w[-1] = 1/3*delta_t
|
||||
return t, w
|
||||
return t, w
|
||||
|
||||
def get_four_point_weights_times(t_max, num_grid_points):
|
||||
if num_grid_points % 4 != 1:
|
||||
raise RuntimeError("four point weights needs grid points ng such that ng = 4*k+1, but git ng={}".format(num_grid_points))
|
||||
# generate mid points
|
||||
t, delta_t = np.linspace(0, t_max, num_grid_points, retstep = True)
|
||||
# equal weights for grid points
|
||||
w = np.empty(num_grid_points, dtype=np.float64)
|
||||
w[0::4] = 2 * 7 * 4 / 90 * delta_t
|
||||
w[1::4] = 1 * 32 * 4 / 90 * delta_t
|
||||
w[2::4] = 1 * 12 * 4 / 90 * delta_t
|
||||
w[3::4] = 1 * 32 * 4 / 90 * delta_t
|
||||
w[0] = 1 * 7 * 4 / 90 * delta_t
|
||||
w[-1] = 1 * 7 * 4 / 90 * delta_t
|
||||
return t, w
|
||||
|
||||
def get_gauss_legendre_weights_times(t_max, num_grid_points):
|
||||
return gquad.gauss_nodes_weights_legendre(n = num_grid_points, low=0, high=t_max)
|
||||
|
||||
def get_sinh_tanh_weights_times(t_max, num_grid_points):
|
||||
"""
|
||||
inspired by Tanh-Sinh High-Precision Quadrature - David H. Bailey
|
||||
|
||||
"""
|
||||
def get_h_of_N(N):
|
||||
"""returns the stepsize h for sinh tanh quad for a given number of points N
|
||||
such that the smallest weights are about 1e-14
|
||||
"""
|
||||
a = 16.12087683080651
|
||||
b = -2.393599730652087
|
||||
c = 6.536936185577097
|
||||
d = -1.012504470475915
|
||||
if N < 4:
|
||||
raise ValueError("only tested for N >= 4")
|
||||
return a * N ** b + c * N ** d
|
||||
|
||||
h = get_h_of_N(num_grid_points)
|
||||
if num_grid_points % 2 != 1:
|
||||
raise RuntimeError("sinh tanh weights needs grid points ng such that ng = 2*k+1, but git ng={}".format(num_grid_points))
|
||||
kmax = (num_grid_points - 1) / 2
|
||||
k = np.arange(0, kmax+1)
|
||||
w = h * np.pi / 2 * np.cosh(k * h) / np.cosh(np.pi / 2 * np.sinh(k * h)) ** 2
|
||||
w = np.hstack((w[-1:0:-1], w))*t_max/2
|
||||
|
||||
tmp = np.pi/2*np.sinh(h*k)
|
||||
y_plus = 1/ (np.exp(tmp) * np.cosh(tmp))
|
||||
t = np.hstack((y_plus[-1:0:-1], (2-y_plus)))*t_max/2
|
||||
return t, w
|
||||
|
||||
|
||||
_WC_ = 2
|
||||
_S_ = 0.6
|
||||
_GAMMA_S_PLUS_1 = gamma(_S_ + 1)
|
||||
def lac(t):
|
||||
"""lorenzian bath correlation function"""
|
||||
return np.exp(- np.abs(t) - 1j * _WC_ * t)
|
||||
|
||||
def oac(t):
|
||||
"""ohmic bath correlation function"""
|
||||
return (1 + 1j*(t))**(-(_S_+1)) * _GAMMA_S_PLUS_1 / np.pi
|
||||
|
||||
def get_rel_diff(corr, t_max, ng, weights_times, ng_fac):
|
||||
t, w = weights_times(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
ng_fine = ng_fac * (ng - 1) + 1
|
||||
tfine = np.linspace(0, t_max, ng_fine)
|
||||
bcf_n_plus = corr(tfine - tfine[0])
|
||||
alpha_k = np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus))
|
||||
|
||||
u_i_all_t = stocproc_c.eig_func_all_interp(delta_t_fac = ng_fac,
|
||||
time_axis = t,
|
||||
alpha_k = alpha_k,
|
||||
weights = w,
|
||||
eigen_val = _eig_val,
|
||||
eigen_vec = _eig_vec)
|
||||
|
||||
u_i_all_ast_s = np.conj(u_i_all_t) # (N_gp, N_ev)
|
||||
num_ev = len(_eig_val)
|
||||
tmp = _eig_val.reshape(1, num_ev) * u_i_all_t # (N_gp, N_ev)
|
||||
recs_bcf = np.tensordot(tmp, u_i_all_ast_s, axes=([1], [1]))
|
||||
|
||||
refc_bcf = np.empty(shape=(ng_fine, ng_fine), dtype=np.complex128)
|
||||
for i in range(ng_fine):
|
||||
idx = ng_fine - 1 - i
|
||||
refc_bcf[:, i] = alpha_k[idx:idx + ng_fine]
|
||||
|
||||
rd = np.abs(recs_bcf - refc_bcf) / np.abs(refc_bcf)
|
||||
return tfine, rd
|
||||
|
|
@ -40,7 +40,7 @@ cpdef cnp.ndarray[DTYPE_CPLX_t, ndim=1] eig_func_interp(unsigned int
|
|||
|
||||
|
||||
|
||||
cpdef cnp.ndarray[DTYPE_CPLX_t, ndim=2] eig_func_all_interp(unsigned int delta_t_fac,
|
||||
cpdef cnp.ndarray[DTYPE_CPLX_t, ndim=2] eig_func_all_interp(unsigned int delta_t_fac,
|
||||
cnp.ndarray[DTYPE_DBL_t, ndim=1] time_axis,
|
||||
cnp.ndarray[DTYPE_CPLX_t, ndim=1] alpha_k,
|
||||
cnp.ndarray[DTYPE_DBL_t, ndim=1] weights,
|
||||
|
|
|
|||
|
|
@ -90,3 +90,5 @@ def auto_correlation_zero(x, s_0_idx = 0):
|
|||
num_samples = x.shape[0]
|
||||
x_s_0 = x[:,s_0_idx].reshape(num_samples,1)
|
||||
return np.mean(x * np.conj(x_s_0), axis = 0), np.mean(x * x_s_0, axis = 0)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@ import numpy as np
|
|||
import math
|
||||
from scipy.special import gamma as gamma_func
|
||||
import scipy.integrate as sp_int
|
||||
from scipy.linalg import eig
|
||||
try:
|
||||
import matplotlib.pyplot as plt
|
||||
except ImportError:
|
||||
|
|
@ -17,88 +18,547 @@ sys.path.insert(0, str(p.parent.parent))
|
|||
|
||||
import stocproc as sp
|
||||
from stocproc import tools
|
||||
from stocproc import method_kle
|
||||
from stocproc import stocproc_c
|
||||
from scipy.integrate import quad
|
||||
import pickle
|
||||
import logging
|
||||
|
||||
def test_weights(plot=False):
|
||||
"""
|
||||
here we check for the correct scaling of the various integration schemas
|
||||
"""
|
||||
def f(x):
|
||||
return x/(1+x**2)
|
||||
|
||||
tm = 10
|
||||
meth = [method_kle.get_mid_point_weights_times,
|
||||
method_kle.get_simpson_weights_times,
|
||||
method_kle.get_four_point_weights_times,
|
||||
method_kle.get_gauss_legendre_weights_times,
|
||||
method_kle.get_sinh_tanh_weights_times]
|
||||
cols = ['r', 'b', 'g', 'm', 'c']
|
||||
errs = [2e-3, 2e-8, 7e-11, 5e-16, 8e-15]
|
||||
I_exact = np.log(tm**2 + 1)/2
|
||||
|
||||
ng = 401
|
||||
for i, _meth in enumerate(meth):
|
||||
t, w = _meth(t_max=tm, num_grid_points=ng)
|
||||
err = abs(I_exact - np.sum(w * f(t)))
|
||||
print(_meth.__name__, err)
|
||||
assert err < errs[i]
|
||||
|
||||
if plot:
|
||||
for i, _meth in enumerate(meth):
|
||||
xdata = []
|
||||
ydata = []
|
||||
for k in np.logspace(0, 2, 50):
|
||||
ng = 4 * int(k) + 1
|
||||
t, w = _meth(t_max=tm, num_grid_points=ng)
|
||||
I = np.sum(w*f(t))
|
||||
xdata.append(ng)
|
||||
ydata.append(abs(I - I_exact))
|
||||
plt.plot(xdata, ydata, marker='o', color=cols[i], label=_meth.__name__)
|
||||
|
||||
plt.legend(loc = 'lower left')
|
||||
plt.xscale('log')
|
||||
plt.yscale('log')
|
||||
plt.grid()
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
def test_solve_fredholm():
|
||||
"""
|
||||
here we compare the fredholm eigenvalue problem solver
|
||||
(which maps the non hermetian eigenvalue problem to a hermetian one)
|
||||
with a straigt forward non hermetian eigenvalue solver
|
||||
"""
|
||||
t_max = 15
|
||||
corr = method_kle.lac
|
||||
meth = [method_kle.get_mid_point_weights_times,
|
||||
method_kle.get_simpson_weights_times,
|
||||
method_kle.get_four_point_weights_times]
|
||||
ng = 41
|
||||
for _meth in meth:
|
||||
t, w = _meth(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w)
|
||||
eval, evec = eig(r*w.reshape(1,-1))
|
||||
max_imag = np.max(np.abs(np.imag(eval)))
|
||||
assert max_imag < 1e-15 # make sure the evals are real
|
||||
|
||||
eval = np.real(eval)
|
||||
idx_sort = np.argsort(eval)[::-1]
|
||||
eval = eval[idx_sort]
|
||||
evec = evec[:,idx_sort]
|
||||
max_eval_diff = np.max(np.abs(eval - _eig_val))
|
||||
assert max_eval_diff < 1e-13 # compare evals with fredholm evals
|
||||
|
||||
for i in range(ng-5):
|
||||
phase = np.arctan2(np.imag(_eig_vec[0,i]), np.real(_eig_vec[0,i]))
|
||||
_eig_vec[:, i] *= np.exp(-1j * phase)
|
||||
norm = np.sqrt(np.sum(w*np.abs(_eig_vec[:, i]) ** 2))
|
||||
assert abs(1-norm) < 1e-14 # fredholm should return vectors with integral norm 1
|
||||
|
||||
phase = np.arctan2(np.imag(evec[0, i]), np.real(evec[0, i]))
|
||||
evec[:, i] *= np.exp(-1j * phase)
|
||||
norm = np.sqrt(np.sum(w* np.abs(evec[:, i])**2))
|
||||
evec[:, i] /= norm
|
||||
|
||||
diff_evec = np.max(np.abs(_eig_vec[:,i] - evec[:,i]))
|
||||
assert diff_evec < 1e-10, "diff_evec {} = {} {}".format(i, diff_evec, _meth.__name__)
|
||||
|
||||
|
||||
def test_compare_weights_in_solve_fredholm_oac():
|
||||
"""
|
||||
here we try to examine which integration weights perform best in order to
|
||||
calculate the eigenfunctions -> well it seems to depend on the situation
|
||||
|
||||
although simpson and gauss-legendre perform well
|
||||
"""
|
||||
t_max = 15
|
||||
corr = method_kle.oac
|
||||
|
||||
ng_ref = 3501
|
||||
|
||||
_meth_ref = method_kle.get_simpson_weights_times
|
||||
t, w = _meth_ref(t_max, ng_ref)
|
||||
|
||||
try:
|
||||
with open("test_fredholm_interpolation.dump", 'rb') as f:
|
||||
ref_data = pickle.load(f)
|
||||
except FileNotFoundError:
|
||||
ref_data = {}
|
||||
key = (tuple(t), tuple(w), corr.__name__)
|
||||
if key in ref_data:
|
||||
eigval_ref, evec_ref = ref_data[key]
|
||||
else:
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
eigval_ref, evec_ref = method_kle.solve_hom_fredholm(r, w)
|
||||
ref_data[key] = eigval_ref, evec_ref
|
||||
with open("test_fredholm_interpolation.dump", 'wb') as f:
|
||||
pickle.dump(ref_data, f)
|
||||
|
||||
method_kle.align_eig_vec(evec_ref)
|
||||
|
||||
ks = [20, 40, 80, 160]
|
||||
|
||||
fig, ax = plt.subplots(ncols=2, nrows=2, sharex=True, sharey=True, figsize=(16, 12))
|
||||
|
||||
ax = ax.flatten()
|
||||
|
||||
lines = []
|
||||
labels = []
|
||||
|
||||
eigvec_ref = []
|
||||
for i in range(ng_ref):
|
||||
eigvec_ref.append(tools.ComplexInterpolatedUnivariateSpline(t, evec_ref[:, i]))
|
||||
|
||||
meth = [method_kle.get_mid_point_weights_times,
|
||||
method_kle.get_simpson_weights_times,
|
||||
method_kle.get_four_point_weights_times,
|
||||
method_kle.get_gauss_legendre_weights_times,
|
||||
method_kle.get_sinh_tanh_weights_times]
|
||||
cols = ['r', 'b', 'g', 'm', 'c']
|
||||
for j, k in enumerate(ks):
|
||||
axc = ax[j]
|
||||
ng = 4 * k + 1
|
||||
|
||||
for i, _meth in enumerate(meth):
|
||||
t, w = _meth(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
method_kle.align_eig_vec(_eig_vec)
|
||||
|
||||
dx = []
|
||||
dy = []
|
||||
dy2 = []
|
||||
|
||||
for l in range(len(_eig_val)):
|
||||
evr = eigvec_ref[l](t)
|
||||
diff = np.abs(_eig_vec[:, l] - evr)
|
||||
dx.append(l)
|
||||
dy.append(np.max(diff))
|
||||
dy2.append(abs(_eig_val[l] - eigval_ref[l]))
|
||||
|
||||
p, = axc.plot(dx, dy, color=cols[i])
|
||||
axc.plot(dx, dy2, color=cols[i], ls='--')
|
||||
if j == 0:
|
||||
lines.append(p)
|
||||
labels.append(_meth.__name__)
|
||||
|
||||
t, w = method_kle.get_simpson_weights_times(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
method_kle.align_eig_vec(_eig_vec)
|
||||
_eig_vec = _eig_vec[1::2, :]
|
||||
t = t[1::2]
|
||||
|
||||
dx = []
|
||||
dy = []
|
||||
dy2 = []
|
||||
|
||||
for l in range(len(_eig_val)):
|
||||
evr = eigvec_ref[l](t)
|
||||
diff = np.abs(_eig_vec[:, l] - evr)
|
||||
dx.append(l)
|
||||
dy.append(np.max(diff))
|
||||
dy2.append(abs(_eig_val[l] - eigval_ref[l]))
|
||||
|
||||
p, = axc.plot(dx, dy, color='lime')
|
||||
p_lam, = axc.plot(list(range(ng)), eigval_ref[:ng], color='k')
|
||||
if j == 0:
|
||||
lines.append(p)
|
||||
labels.append('simp 2nd')
|
||||
lines.append(p_lam)
|
||||
labels.append('abs(lamnda_i)')
|
||||
|
||||
axc.grid()
|
||||
axc.set_yscale('log')
|
||||
axc.set_xlim([0, 100])
|
||||
axc.set_ylim([1e-15, 10])
|
||||
axc.set_title("ng {}".format(ng))
|
||||
|
||||
fig.suptitle("use ref with ng_ref {} and method '{}'".format(ng_ref, _meth_ref.__name__))
|
||||
fig.legend(handles=lines, labels=labels, ncol=3, loc='lower center')
|
||||
fig.subplots_adjust(bottom=0.15)
|
||||
plt.show()
|
||||
|
||||
def test_compare_weights_in_solve_fredholm_lac():
|
||||
"""
|
||||
here we try to examine which integration weights perform best in order to
|
||||
calculate the eigenfunctions -> well it seems to depend on the situation
|
||||
|
||||
although simpson and gauss-legendre perform well
|
||||
"""
|
||||
t_max = 15
|
||||
corr = method_kle.lac
|
||||
|
||||
ng_ref = 3501
|
||||
|
||||
_meth_ref = method_kle.get_simpson_weights_times
|
||||
t, w = _meth_ref(t_max, ng_ref)
|
||||
|
||||
try:
|
||||
with open("test_fredholm_interpolation.dump", 'rb') as f:
|
||||
ref_data = pickle.load(f)
|
||||
except FileNotFoundError:
|
||||
ref_data = {}
|
||||
key = (tuple(t), tuple(w), corr.__name__)
|
||||
if key in ref_data:
|
||||
eigval_ref, evec_ref = ref_data[key]
|
||||
else:
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
eigval_ref, evec_ref = method_kle.solve_hom_fredholm(r, w)
|
||||
ref_data[key] = eigval_ref, evec_ref
|
||||
with open("test_fredholm_interpolation.dump", 'wb') as f:
|
||||
pickle.dump(ref_data, f)
|
||||
|
||||
|
||||
method_kle.align_eig_vec(evec_ref)
|
||||
|
||||
ks = [20,40,80,160]
|
||||
|
||||
fig, ax = plt.subplots(ncols=2, nrows=2, sharex=True, sharey=True, figsize=(16,12))
|
||||
|
||||
ax = ax.flatten()
|
||||
|
||||
lines = []
|
||||
labels = []
|
||||
|
||||
eigvec_ref = []
|
||||
for i in range(ng_ref):
|
||||
eigvec_ref.append(tools.ComplexInterpolatedUnivariateSpline(t, evec_ref[:, i]))
|
||||
|
||||
meth = [method_kle.get_mid_point_weights_times,
|
||||
method_kle.get_simpson_weights_times,
|
||||
method_kle.get_four_point_weights_times,
|
||||
method_kle.get_gauss_legendre_weights_times,
|
||||
method_kle.get_sinh_tanh_weights_times]
|
||||
cols = ['r', 'b', 'g', 'm', 'c']
|
||||
for j, k in enumerate(ks):
|
||||
axc = ax[j]
|
||||
ng = 4*k+1
|
||||
|
||||
for i, _meth in enumerate(meth):
|
||||
t, w = _meth(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
method_kle.align_eig_vec(_eig_vec)
|
||||
|
||||
dx = []
|
||||
dy = []
|
||||
dy2 = []
|
||||
|
||||
for l in range(len(_eig_val)):
|
||||
evr = eigvec_ref[l](t)
|
||||
diff = np.abs(_eig_vec[:,l] - evr)
|
||||
dx.append(l)
|
||||
dy.append(np.max(diff))
|
||||
dy2.append(abs(_eig_val[l] - eigval_ref[l]))
|
||||
|
||||
p, = axc.plot(dx, dy, color=cols[i])
|
||||
axc.plot(dx, dy2, color=cols[i], ls='--')
|
||||
if j == 0:
|
||||
lines.append(p)
|
||||
labels.append(_meth.__name__)
|
||||
|
||||
|
||||
t, w = method_kle.get_simpson_weights_times(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
method_kle.align_eig_vec(_eig_vec)
|
||||
_eig_vec = _eig_vec[1::2, :]
|
||||
t = t[1::2]
|
||||
|
||||
dx = []
|
||||
dy = []
|
||||
dy2 = []
|
||||
|
||||
for l in range(len(_eig_val)):
|
||||
evr = eigvec_ref[l](t)
|
||||
diff = np.abs(_eig_vec[:, l] - evr)
|
||||
dx.append(l)
|
||||
dy.append(np.max(diff))
|
||||
dy2.append(abs(_eig_val[l] - eigval_ref[l]))
|
||||
|
||||
p, = axc.plot(dx, dy, color='lime')
|
||||
p_lam, = axc.plot(list(range(ng)), eigval_ref[:ng], color='k')
|
||||
if j == 0:
|
||||
lines.append(p)
|
||||
labels.append('simp 2nd')
|
||||
lines.append(p_lam)
|
||||
labels.append('abs(lamnda_i)')
|
||||
|
||||
axc.grid()
|
||||
axc.set_yscale('log')
|
||||
axc.set_xlim([0,100])
|
||||
axc.set_ylim([1e-5, 10])
|
||||
axc.set_title("ng {}".format(ng))
|
||||
|
||||
|
||||
fig.suptitle("use ref with ng_ref {} and method '{}'".format(ng_ref, _meth_ref.__name__))
|
||||
fig.legend(handles=lines, labels=labels, ncol=3, loc='lower center')
|
||||
fig.subplots_adjust(bottom=0.15)
|
||||
plt.show()
|
||||
|
||||
def test_fredholm_eigvec_interpolation():
|
||||
"""
|
||||
for ohmic sd : use 4 point and integral interpolation
|
||||
for lorentzian : use simpson and spline interpolation
|
||||
"""
|
||||
t_max = 15
|
||||
corr = method_kle.lac
|
||||
#corr = method_kle.oac
|
||||
|
||||
ng_ref = 3501
|
||||
|
||||
_meth_ref = method_kle.get_simpson_weights_times
|
||||
t, w = _meth_ref(t_max, ng_ref)
|
||||
|
||||
try:
|
||||
with open("test_fredholm_interpolation.dump", 'rb') as f:
|
||||
ref_data = pickle.load(f)
|
||||
except FileNotFoundError:
|
||||
ref_data = {}
|
||||
key = (tuple(t), tuple(w), corr.__name__)
|
||||
if key in ref_data:
|
||||
eigval_ref, evec_ref = ref_data[key]
|
||||
else:
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
eigval_ref, evec_ref = method_kle.solve_hom_fredholm(r, w)
|
||||
ref_data[key] = eigval_ref, evec_ref
|
||||
with open("test_fredholm_interpolation.dump", 'wb') as f:
|
||||
pickle.dump(ref_data, f)
|
||||
|
||||
method_kle.align_eig_vec(evec_ref)
|
||||
t_ref = t
|
||||
|
||||
eigvec_ref = []
|
||||
for l in range(ng_ref):
|
||||
eigvec_ref.append(tools.ComplexInterpolatedUnivariateSpline(t, evec_ref[:, l]))
|
||||
|
||||
meth = [method_kle.get_mid_point_weights_times,
|
||||
method_kle.get_simpson_weights_times,
|
||||
method_kle.get_four_point_weights_times,
|
||||
method_kle.get_gauss_legendre_weights_times]
|
||||
cols = ['r', 'b', 'g', 'm', 'c']
|
||||
|
||||
def my_intp(ti, corr, w, t, u, lam):
|
||||
return np.sum(corr(ti - t) * w * u) / lam
|
||||
|
||||
fig, ax = plt.subplots(ncols=2, nrows=2, sharex=True, sharey=True, figsize=(16,12))
|
||||
ax = ax.flatten()
|
||||
|
||||
ks = [20,40,80,160]
|
||||
|
||||
lns, lbs = [], []
|
||||
|
||||
for ik, k in enumerate(ks):
|
||||
axc = ax[ik]
|
||||
|
||||
ng = 4*k+1
|
||||
for i, _meth in enumerate(meth):
|
||||
print(ik, i)
|
||||
t, w = _meth(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w)
|
||||
method_kle.align_eig_vec(_eig_vec)
|
||||
|
||||
eigvec_intp = []
|
||||
for l in range(ng):
|
||||
eigvec_intp.append(tools.ComplexInterpolatedUnivariateSpline(t, _eig_vec[:, l]))
|
||||
|
||||
ydata_fixed = []
|
||||
ydata_spline = []
|
||||
ydata_integr_intp = []
|
||||
xdata = np.arange(min(ng, 100))
|
||||
|
||||
for idx in xdata:
|
||||
evr = eigvec_ref[idx](t)
|
||||
ydata_fixed.append(np.max(np.abs(_eig_vec[:,idx] - evr)))
|
||||
ydata_spline.append(np.max(np.abs(eigvec_intp[idx](t_ref) - evec_ref[:,idx])))
|
||||
uip = np.asarray([my_intp(ti, corr, w, t, _eig_vec[:,idx], _eig_val[idx]) for ti in t_ref])
|
||||
ydata_integr_intp.append(np.max(np.abs(uip - evec_ref[:,idx])))
|
||||
|
||||
p1, = axc.plot(xdata, ydata_fixed, color=cols[i], label=_meth.__name__)
|
||||
p2, = axc.plot(xdata, ydata_spline, color=cols[i], ls='--')
|
||||
p3, = axc.plot(xdata, ydata_integr_intp, color=cols[i], alpha = 0.5)
|
||||
if ik == 0:
|
||||
lns.append(p1)
|
||||
lbs.append(_meth.__name__)
|
||||
|
||||
if ik == 0:
|
||||
lines = [p1,p2,p3]
|
||||
labels = ['fixed', 'spline', 'integral interp']
|
||||
axc.set_yscale('log')
|
||||
axc.set_title("ng {}".format(ng))
|
||||
axc.set_xlim([0,100])
|
||||
axc.grid()
|
||||
|
||||
|
||||
fig.legend(lines, labels, loc = "lower right", ncol=3)
|
||||
fig.legend(lns, lbs, loc="lower left", ncol=2)
|
||||
plt.subplots_adjust(bottom = 0.15)
|
||||
|
||||
plt.savefig("test_fredholm_eigvec_interpolation_{}_.pdf".format(corr.__name__))
|
||||
plt.show()
|
||||
|
||||
def test_oac_interpolation():
|
||||
"""
|
||||
"""
|
||||
t_max = 15
|
||||
corr = method_kle.oac
|
||||
|
||||
# ng_ref = 3501
|
||||
# _meth_ref = method_kle.get_simpson_weights_times
|
||||
# t, w = _meth_ref(t_max, ng_ref)
|
||||
#
|
||||
# try:
|
||||
# with open("test_fredholm_interpolation.dump", 'rb') as f:
|
||||
# ref_data = pickle.load(f)
|
||||
# except FileNotFoundError:
|
||||
# ref_data = {}
|
||||
# key = (tuple(t), tuple(w), corr.__name__)
|
||||
# if key in ref_data:
|
||||
# eigval_ref, evec_ref = ref_data[key]
|
||||
# else:
|
||||
# r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
# eigval_ref, evec_ref = method_kle.solve_hom_fredholm(r, w)
|
||||
# ref_data[key] = eigval_ref, evec_ref
|
||||
# with open("test_fredholm_interpolation.dump", 'wb') as f:
|
||||
# pickle.dump(ref_data, f)
|
||||
#
|
||||
# method_kle.align_eig_vec(evec_ref)
|
||||
# t_ref = t
|
||||
#
|
||||
# eigvec_ref = []
|
||||
# for l in range(ng_ref):
|
||||
# eigvec_ref.append(tools.ComplexInterpolatedUnivariateSpline(t, evec_ref[:, l]))
|
||||
#
|
||||
meth = method_kle.get_four_point_weights_times
|
||||
def my_intp(ti, corr, w, t, u, lam):
|
||||
return np.sum(u * corr(ti - t) * w) / lam
|
||||
|
||||
k = 40
|
||||
ng = 4*k+1
|
||||
|
||||
t, w = meth(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w)
|
||||
method_kle.align_eig_vec(_eig_vec)
|
||||
|
||||
ngfac = 4
|
||||
|
||||
tfine = np.linspace(0, t_max, (ng-1)*ngfac+1)
|
||||
|
||||
#ui_fine = np.empty(shape=(ng//2, len(tfine)), dtype = np.complex128)
|
||||
|
||||
i = 40
|
||||
|
||||
ui_fine = np.asarray([my_intp(ti, corr, w, t, _eig_vec[:,i], _eig_val[i]) for ti in tfine])
|
||||
|
||||
bcf_n_plus = corr(tfine - tfine[0])
|
||||
alpha_k = np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus))
|
||||
|
||||
# ui_all_t = stocproc_c.eig_func_all_interp(delta_t_fac=ngfac,
|
||||
# time_axis=t,
|
||||
# alpha_k=alpha_k,
|
||||
# weights=w,
|
||||
# eigen_val=_eig_val[:ng//2],
|
||||
# eigen_vec=_eig_vec[:, :ng//2])
|
||||
ui_fine2 = stocproc_c.eig_func_interp(delta_t_fac = ngfac,
|
||||
time_axis = t,
|
||||
alpha_k = alpha_k,
|
||||
weights = w,
|
||||
eigen_val = _eig_val[i],
|
||||
eigen_vec = _eig_vec[:,i])
|
||||
plt.plot(tfine, np.abs(ui_fine - ui_fine2)*_eig_val[i])
|
||||
plt.yscale('log')
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
def test_opt_fredh():
|
||||
t_max = 15
|
||||
corr = method_kle.lac
|
||||
|
||||
ng = 1001
|
||||
|
||||
t, w = method_kle.get_simpson_weights_times(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
eigval, eigvec = method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
idx = 0
|
||||
method_kle.opt_fredholm(corr, eigval[idx], t, eigvec[:,idx], method_kle.get_simpson_weights_times, ntimes=3)
|
||||
|
||||
|
||||
def test_solve_fredholm2():
|
||||
_WC_ = 2
|
||||
def lac(t):
|
||||
return np.exp(- np.abs(t) - 1j*_WC_*t)
|
||||
|
||||
t_max = 10
|
||||
ng_fac = 1
|
||||
idx_ = 102
|
||||
|
||||
for ng in range(11, 450, 30):
|
||||
t, w = sp.method_kle.get_mid_point_weights_times(t_max, ng)
|
||||
r = lac(t.reshape(-1,1)-t.reshape(1,-1))
|
||||
_eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
corr = method_kle.lac
|
||||
|
||||
ng_fac = 4
|
||||
ng_fine = ng_fac * (ng - 1) + 1
|
||||
tfine = np.linspace(0, t_max, ng_fine)
|
||||
bcf_n_plus = lac(tfine - tfine[0])
|
||||
alpha_k = np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus))
|
||||
for ng in [51]:
|
||||
tfine, rd = method_kle.get_rel_diff(corr = corr, t_max=t_max, ng=ng,
|
||||
weights_times=method_kle.get_mid_point_weights_times, ng_fac=ng_fac)
|
||||
print(np.max(rd))
|
||||
plt.imshow(np.log(rd))
|
||||
#plt.plot(tfine, rd[:, idx_])
|
||||
|
||||
u_i_all_t = stocproc_c.eig_func_all_interp(delta_t_fac=ng_fac,
|
||||
time_axis=t,
|
||||
alpha_k=alpha_k,
|
||||
weights=w,
|
||||
eigen_val=_eig_val,
|
||||
eigen_vec=_eig_vec)
|
||||
|
||||
u_i_all_ast_s = np.conj(u_i_all_t) # (N_gp, N_ev)
|
||||
num_ev = len(_eig_val)
|
||||
tmp = _eig_val.reshape(1, num_ev) * u_i_all_t # (N_gp, N_ev)
|
||||
recs_bcf = np.tensordot(tmp, u_i_all_ast_s, axes=([1], [1]))
|
||||
|
||||
refc_bcf = np.empty(shape=(ng_fine, ng_fine), dtype=np.complex128)
|
||||
for i in range(ng_fine):
|
||||
idx = ng_fine - 1 - i
|
||||
refc_bcf[:, i] = alpha_k[idx:idx + ng_fine]
|
||||
|
||||
err = np.max(np.abs(recs_bcf - refc_bcf) / np.abs(refc_bcf))
|
||||
plt.plot(ng, err, marker='o', color='r')
|
||||
|
||||
ng += 2
|
||||
t, w = sp.method_kle.get_simpson_weights_times(t_max, ng)
|
||||
r = lac(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
_eig_val2, _eig_vec2 = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0)
|
||||
#print(np.max(np.abs(_eig_val - _eig_val2)), np.max(np.abs(_eig_vec - _eig_vec2)))
|
||||
_eig_val, _eig_vec = _eig_val2, _eig_vec2
|
||||
|
||||
ng_fac = 4
|
||||
ng_fine = ng_fac * (ng - 1) + 1
|
||||
tfine = np.linspace(0, t_max, ng_fine)
|
||||
bcf_n_plus = lac(tfine - tfine[0])
|
||||
alpha_k = np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus))
|
||||
|
||||
u_i_all_t2 = stocproc_c.eig_func_all_interp(delta_t_fac=ng_fac,
|
||||
time_axis=t,
|
||||
alpha_k=alpha_k,
|
||||
weights=w,
|
||||
eigen_val=_eig_val,
|
||||
eigen_vec=_eig_vec)
|
||||
#print(np.max(np.abs(u_i_all_t - u_i_all_t2)))
|
||||
u_i_all_t = u_i_all_t2
|
||||
#print()
|
||||
#tfine, rd = method_kle.get_rel_diff(corr=corr, t_max=t_max, ng=ng,
|
||||
# weights_times=method_kle.get_simpson_weights_times, ng_fac=ng_fac)
|
||||
#plt.plot(tfine, rd[:, idx_])
|
||||
|
||||
|
||||
u_i_all_ast_s = np.conj(u_i_all_t) # (N_gp, N_ev)
|
||||
num_ev = len(_eig_val)
|
||||
tmp = _eig_val.reshape(1, num_ev) * u_i_all_t # (N_gp, N_ev)
|
||||
recs_bcf = np.tensordot(tmp, u_i_all_ast_s, axes=([1], [1]))
|
||||
|
||||
refc_bcf = np.empty(shape=(ng_fine, ng_fine), dtype=np.complex128)
|
||||
for i in range(ng_fine):
|
||||
idx = ng_fine - 1 - i
|
||||
refc_bcf[:, i] = alpha_k[idx:idx + ng_fine]
|
||||
|
||||
err2 = np.max(np.abs(recs_bcf - refc_bcf) / np.abs(refc_bcf))
|
||||
plt.plot(ng, err2, marker='o', color='k')
|
||||
print(err, err2)
|
||||
|
||||
plt.yscale('log')
|
||||
plt.grid()
|
||||
#plt.yscale('log')
|
||||
#plt.grid()
|
||||
plt.show()
|
||||
|
||||
def test_solve_fredholm_reconstr_ac():
|
||||
|
|
@ -191,6 +651,14 @@ def test_solve_fredholm_interp_eigenfunc(plot=False):
|
|||
|
||||
|
||||
if __name__ == "__main__":
|
||||
# test_weights(plot=True)
|
||||
# test_solve_fredholm()
|
||||
# test_compare_weights_in_solve_fredholm_oac()
|
||||
# test_compare_weights_in_solve_fredholm_lac()
|
||||
# test_fredholm_eigvec_interpolation()
|
||||
test_oac_interpolation()
|
||||
# test_opt_fredh()
|
||||
# test_solve_fredholm()
|
||||
# test_solve_fredholm_reconstr_ac()
|
||||
test_solve_fredholm_interp_eigenfunc(plot=True)
|
||||
# test_solve_fredholm_interp_eigenfunc(plot=True)
|
||||
pass
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue