hiro/stocproc
1
0
Fork 0
mirror of https://github.com/vale981/stocproc synced 2025-03-06 02:01:41 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

merged selected stuff from homework

Browse source
This commit is contained in:
Richard Hartmann 2016-10-28 10:12:32 +02:00
parent 6b288ee48b
commit 263a69cfc7
3 changed files with 24 additions and 16 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
stocproc/class_stocproc_kle.py
View file

@ -598,11 +598,11 @@ def recons_corr_and_get_bcf(T, ng, w, eig_val, eig_vec, bcf):
delta_t_fac = 2 delta_t_fac = 2
N1 = ng N1 = ng
N2 = delta_t_fac * (N1 - 1) + 1 N2 = delta_t_fac * (N1 - 1) + 1
alpha_k = bcf(np.linspace(-T, T, 2*N2 - 1)) bcf_n_plus = bcf(np.linspace(0, T, N2))
bcf_n = 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 = delta_t_fac, u_i_all_t = stocproc_c.eig_func_all_interp(delta_t_fac = delta_t_fac,
time_axis = np.linspace(0, T, N1), time_axis = np.linspace(0, T, N1),
alpha_k = alpha_k, alpha_k = bcf_n,
weights = w, weights = w,
eigen_val = eig_val, eigen_val = eig_val,
eigen_vec = eig_vec) eigen_vec = eig_vec)
@ -615,12 +615,12 @@ def recons_corr_and_get_bcf(T, ng, w, eig_val, eig_vec, bcf):
refc_bcf = np.empty(shape=(N2,N2), dtype = np.complex128) refc_bcf = np.empty(shape=(N2,N2), dtype = np.complex128)
for i in range(N2): for i in range(N2):
idx = N2-1-i idx = N2-1-i
refc_bcf[:,i] = alpha_k[idx:idx+N2] refc_bcf[:,i] = bcf_n[idx:idx+N2]
return recs_bcf, refc_bcf return recs_bcf, refc_bcf
def auto_grid_points(r_tau, t_max, tol = 1e-8, err_method = max_error, name = 'mid_point', sig_min = 1e-4, verbose=1): def auto_grid_points(r_tau, t_max, tol = 1e-3, err_method = max_rel_error, name = 'simpson', sig_min = 1e-6, verbose=1):
err = 1 err = 1
c = 2 c = 2
seed = None seed = None
@ -632,7 +632,6 @@ def auto_grid_points(r_tau, t_max, tol = 1e-8, err_method = max_error, name = 'm
c *= 2 c *= 2
ng = 2*c + 1 ng = 2*c + 1
ng_fine = ng*2-1 ng_fine = ng*2-1
t_fine = np.linspace(0, t_max, ng_fine)
if verbose == 1: if verbose == 1:
print("ng:{} new proc ({}) ... ".format(ng, name), end='') print("ng:{} new proc ({}) ... ".format(ng, name), end='')
sys.stdout.flush() sys.stdout.flush()
@ -665,23 +664,27 @@ def auto_grid_points(r_tau, t_max, tol = 1e-8, err_method = max_error, name = 'm
print("c_high", c_high) print("c_high", c_high)
c = (c_low + c_high) // 2 c = (c_low + c_high) // 2
ng = 2*c + 1 ng = 2*c + 1
ng_fine = ng * 2 - 1
if verbose > 1: if verbose > 1:
print("c", c)
print("ng", ng) print("ng", ng)
ng_fine = ng*2-1 print("ng_fine", ng_fine)
t_fine = np.linspace(0, t_max, ng_fine)
if verbose == 1: if verbose == 1:
print("ng:{} new proc ({}) ... ".format(ng, name), end='') print("ng:{} new proc ({}) ... ".format(ng, name), end='')
sys.stdout.flush() sys.stdout.flush()
if verbose > 1: if verbose > 1:
print("new process with {} weights ...".format(name)) print("new process with {} weights ({} points)...".format(name, ng))
stoc_proc = StocProc.new_instance_by_name(name, r_tau, t_max, ng, seed, sig_min, verbose-1) stoc_proc = StocProc.new_instance_by_name(name, r_tau, t_max, ng, seed, sig_min, verbose-1)
if verbose > 1: if verbose > 1:
print("reconstruct correlation function ({} points)...".format(ng_fine)) print("reconstruct correlation function ({} points)...".format(ng_fine))
r_t_s = stoc_proc.recons_corr(t_fine) r_t_s, r_t_s_exact = recons_corr_and_get_bcf(T = t_max,
if verbose > 1: ng = ng,
print("calculate exact correlation function ...") w = stoc_proc._w,
r_t_s_exact = r_tau(t_fine.reshape(ng_fine,1) - t_fine.reshape(1, ng_fine)) eig_val = stoc_proc._eig_val,
eig_vec = stoc_proc._eig_vec,
bcf = r_tau)
if verbose > 1: if verbose > 1:
print("calculate error using {} ...".format(err_method_name)) print("calculate error using {} ...".format(err_method_name))
err = np.max(err_method(r_t_s, r_t_s_exact)) err = np.max(err_method(r_t_s, r_t_s_exact))

2
stocproc/method_fft.py
View file

@ -185,7 +185,7 @@ def get_dt_for_accurate_interpolation(t_max, tol, ft_ref):
return t_max/(N/sub_sampl) return t_max/(N/sub_sampl)
N*=2 N*=2
def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, tmax, a, b, ft_ref, N_max = 2**15, method='simps'): def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, tmax, a, b, ft_ref, N_max = 2**20, method='simps'):
N = get_N_for_accurate_fourier_integral(integrand, a, b, N = get_N_for_accurate_fourier_integral(integrand, a, b,
t_max = tmax, t_max = tmax,
tol = intgr_tol, tol = intgr_tol,

7
tests/test_method_fft.py
View file

@ -14,6 +14,11 @@ sys.path.insert(0, str(p.parent.parent))
import stocproc as sp import stocproc as sp
import logging
logging.basicConfig(level=logging.DEBUG)
def test_find_integral_boundary(): def test_find_integral_boundary():
def f(x): def f(x):
return np.exp(-(x)**2) return np.exp(-(x)**2)
@ -286,7 +291,7 @@ def test_get_N_for_accurate_fourier_integral():
bcf_ref = lambda t: gamma_func(s + 1) * wc**(s+1) * (1 + 1j*wc * t)**(-(s+1)) bcf_ref = lambda t: gamma_func(s + 1) * wc**(s+1) * (1 + 1j*wc * t)**(-(s+1))
a,b = sp.method_fft.find_integral_boundary_auto(integrand=intg, tol=1e-12, ref_val=1) a,b = sp.method_fft.find_integral_boundary_auto(integrand=intg, tol=1e-12, ref_val=1)
N = sp.method_fft.get_N_for_accurate_fourier_integral(intg, a, b, t_max=40, tol=1e-3, ft_ref=bcf_ref, N_max = 2**15, method='simps') N = sp.method_fft.get_N_for_accurate_fourier_integral(intg, a, b, t_max=40, tol=1e-3, ft_ref=bcf_ref, N_max = 2**20, method='simps')
print(N) print(N)
def test_get_dt_for_accurate_interpolation(): def test_get_dt_for_accurate_interpolation():