From 56bcee01d74f6d7687d3ca1e14fd9c1402f33b9b Mon Sep 17 00:00:00 2001 From: Richard Hartmann Date: Wed, 10 May 2017 09:51:01 +0200 Subject: [PATCH] SP TanhSinh works --- stocproc/__init__.py | 1 + stocproc/method_fft.py | 2 +- stocproc/stocproc.py | 63 ++++++++++++++++++++++------------------ tests/test_method_fft.py | 42 +++++++++++++++++++++------ 4 files changed, 69 insertions(+), 39 deletions(-) diff --git a/stocproc/__init__.py b/stocproc/__init__.py index d19f663..cd63846 100644 --- a/stocproc/__init__.py +++ b/stocproc/__init__.py @@ -59,5 +59,6 @@ if sys.version_info.major < 3: try: from .stocproc import StocProc_FFT from .stocproc import StocProc_KLE + from .stocproc import StocProc_TanhSinh except ImportError: pass diff --git a/stocproc/method_fft.py b/stocproc/method_fft.py index 3e8797b..de6a7cf 100644 --- a/stocproc/method_fft.py +++ b/stocproc/method_fft.py @@ -170,7 +170,7 @@ def fourier_integral_TanhSinh_with_feedback(integrand, w_max, tau, h, kmax): I = [] feed_back = "ok" for ti in tau: - r = weights_k[0] * integrand(tmp1)*np.exp(-1j*tmp1*ti) + r = weights_k[0] * integrand(tmp1) * np.exp(-1j * tmp1 * ti) for i in range(1, kmax+1): if (y_k[i] * tmp1) == 0: log.debug("y_k is 0") diff --git a/stocproc/stocproc.py b/stocproc/stocproc.py index 6588812..680e5f4 100644 --- a/stocproc/stocproc.py +++ b/stocproc/stocproc.py @@ -462,7 +462,7 @@ class StocProc_TanhSinh(_absStocProc): h = 0.1 k = 15 kstep = 5 - conv_fac = 2 + conv_fac = 1 old_d = None log.info("find h and kmax for TanhSinh integration ...") while True: @@ -476,7 +476,7 @@ class StocProc_TanhSinh(_absStocProc): if fb == 'ok': k += kstep else: - log.info("lowest diff with h {:.3e}: {:.3e} < tol ({:.3e}) -> new h {:.3e}".format(h, d, intgr_tol, h/2)) + log.info("lowest diff with h {:.3e}: {:.3e} < tol ({:.3e}) -> new h {:.3e}".format(h, d[0], intgr_tol, h/2)) h /= 2 k = 15 kstep *= 2 @@ -511,41 +511,46 @@ class StocProc_TanhSinh(_absStocProc): wk = [method_fft.wk(h, ki) for ki in range(1, k+1)] wk = np.hstack((wk[::-1], [method_fft.wk(h, 0)], wk)) - yk = [method_fft.yk(h, ki) for ki in range(1, k+1)] + yk = np.asarray([method_fft.yk(h, ki) for ki in range(1, k+1)]) + tmp1 = wmax/2 + self.omega_k = np.hstack( (yk[::-1] * tmp1, tmp1, (2 - yk) * tmp1)) + self.fl = np.sqrt(tmp1*wk*spectral_density(self.omega_k)/np.pi) - tmp1 = wmax / 2 - #sd = np.hstack((spectral_density(yk * tmp1, )) + super().__init__(t_max=t_max, + num_grid_points=N, + seed=seed, + scale=scale) + @staticmethod + def get_key(t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2): + return bcf_ref, t_max, intgr_tol, intpl_tol - self.yl = wk*spectral_density() + def __getstate__(self): + return self.fl, self.omega_k, self.num_grid_points, self.t_max, self._seed, self.scale, self.key + def __setstate__(self, state): + self.fl, self.omega_k, num_grid_points, t_max, seed, scale, self.key = state + super().__init__(t_max=t_max, + num_grid_points=num_grid_points, + seed=seed, + scale=scale) + def calc_z(self, y): + r"""calculate + + .. math:: + Z(t_l) = sum_k \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k e^{-\i \omega_k t_l} + """ + z = np.empty(shape=self.num_grid_points, dtype=np.complex128) + for i, ti in enumerate(self.t): + z[i] = np.sum(self.fl*y*np.exp(-1j*self.omega_k*ti)) + + return z - - # assert abs(2 * np.pi - N * dx * dt) < 1e-12 - # - # print("Fourier Integral Boundaries: [{:.3e}, {:.3e}]".format(a, b)) - # print("Number of Nodes : {}".format(N)) - # print("yields dx : {:.3e}".format(dx)) - # print("yields dt : {:.3e}".format(dt)) - # print("yields t_max : {:.3e}".format((N - 1) * dt)) - # - # num_grid_points = int(np.ceil(t_max / dt)) + 1 - # - # assert num_grid_points <= N - # - # t_max = (num_grid_points - 1) * dt - # - # super().__init__(t_max=t_max, - # num_grid_points=num_grid_points, - # seed=seed, - # scale=scale) - # - # self.yl = spectral_density(dx * np.arange(N) + a + dx / 2) * dx / np.pi - # self.yl = np.sqrt(self.yl) - # self.omega_min_correction = np.exp(-1j * (a + dx / 2) * self.t) # self.t is from the parent class \ No newline at end of file + def get_num_y(self): + return len(self.fl) \ No newline at end of file diff --git a/tests/test_method_fft.py b/tests/test_method_fft.py index cd4da69..f628c23 100644 --- a/tests/test_method_fft.py +++ b/tests/test_method_fft.py @@ -378,17 +378,41 @@ def test_calc_abN(): +def test_SP_TanhSinh(): + s = -0.5 + wc = 5 + tmax = 25 + + sd = lambda w: w**s*np.exp(-w/wc) + bcf = lambda tau: (wc/(1+1j*wc*tau))**(s+1)*gamma_func(s+1)/np.pi + + _sp = sp.StocProc_TanhSinh(spectral_density=sd, t_max=tmax, bcf_ref=bcf, intgr_tol=1e-2, intpl_tol=1e-2, seed=0) + + d_tol = [0.09, 0.05] + for j, N in enumerate([1000, 5000]): + t = np.linspace(0, tmax, 500) + idx = 200 + c = 0 + for i in range(N): + _sp.new_process() + zt = _sp(t) + c += zt*np.conj(zt[idx]) + c /= N + + d = np.max(np.abs(c - bcf(t-t[idx]))) + print(d) + assert d < d_tol[j] - if __name__ == "__main__": logging.basicConfig(level=logging.INFO) #logging.basicConfig(level=logging.DEBUG) - test_find_integral_boundary() - test_fourier_integral_finite_boundary() - test_fourier_integral_infinite_boundary(plot=False) - test_get_N_a_b_for_accurate_fourier_integral() - test_get_N_a_b_for_accurate_fourier_integral_b_only() - test_get_dt_for_accurate_interpolation() - test_sclicing() - test_calc_abN() \ No newline at end of file + # test_find_integral_boundary() + # test_fourier_integral_finite_boundary() + # test_fourier_integral_infinite_boundary(plot=False) + # test_get_N_a_b_for_accurate_fourier_integral() + # test_get_N_a_b_for_accurate_fourier_integral_b_only() + # test_get_dt_for_accurate_interpolation() + # test_sclicing() + # test_calc_abN() + test_SP_TanhSinh() \ No newline at end of file