From 65f4fb4eeef29dd608d30c92302d2f498397031d Mon Sep 17 00:00:00 2001 From: Richard Hartmann Date: Tue, 9 May 2017 18:21:34 +0200 Subject: [PATCH] work in progress on TanhSinh SP for singular SD --- stocproc/method_fft.py | 98 +++++++++++++++++++------ stocproc/stocproc.py | 163 +++++++++++++++++++++++++++++++++-------- 2 files changed, 207 insertions(+), 54 deletions(-) diff --git a/stocproc/method_fft.py b/stocproc/method_fft.py index de741e4..3e8797b 100644 --- a/stocproc/method_fft.py +++ b/stocproc/method_fft.py @@ -10,6 +10,7 @@ from __future__ import division, print_function #from functools import lru_cache import fcSpline import logging +import mpmath import numpy as np from numpy.fft import rfft as np_rfft from scipy.integrate import quad @@ -131,6 +132,64 @@ def fourier_integral_midpoint(integrand, a, b, N): #log.debug("yields d_x={:.3e}, d_k={:.3e} kmax={:.3e}".format(delta_x, delta_k, tau[-1])) return tau, delta_x*np.exp(-1j*tau*(a+delta_x/2))*fft_vals +def wk(h, k): + return float(0.5*mpmath.pi*h*mpmath.cosh(k*h)/(mpmath.cosh(0.5*mpmath.pi*mpmath.sinh(k*h))**2)) + +def yk(h, k): + return float(1 / (mpmath.exp(mpmath.pi / 2 * mpmath.sinh(k * h)) * mpmath.cosh(mpmath.pi / 2 * mpmath.sinh(k * h)))) + +def fourier_integral_TanhSinh(integrand, w_max, tau, h, kmax): + I, feed_back = fourier_integral_TanhSinh_with_feedback(integrand, w_max, tau, h, kmax) + return I + +def fourier_integral_TanhSinh_with_feedback(integrand, w_max, tau, h, kmax): + """ + + integrate from [0, w_max] the function + integrand*exp(-1j*w*ti) for ti = dt*n, n in [0, N] + + w = w_max/2 (x + 1) # maps the integral from [-1,1] to the integral [a, b] + + weights_k = (0.5 pi cosh(kh)/(cosh(0.5 pi sinh(kh))**2) = weights_minus_k + x_k = 1-y_k = -x_minus_k + y_k = 1/( exp(pi/2 sinh(kh)) cosh(pi/2 np.sinh(kh))) + + I = sum_k weights_k * (b-a)/2 * (integrand(w(x_k)) + integrand(w(x_minus_k))) + + :param integrand: + :param a: + :param b: + :param dt: + :param N: + :return: + """ + k_list = np.arange(kmax+1) + weights_k = [wk(h, ki) for ki in k_list] + y_k = [yk(h, ki) for ki in k_list] + tmp1 = w_max/2 + I = [] + feed_back = "ok" + for ti in tau: + 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") + feed_back = "max kmax reached" + break + + r_tmp = weights_k[i] * ( integrand(y_k[i] * tmp1) * np.exp(-1j*y_k[i] * tmp1*ti) + + integrand((2-y_k[i]) * tmp1) * np.exp(-1j*(2-y_k[i]) * tmp1*ti)) + if np.isnan(r_tmp): + log.debug("integrand yields nan at {} or {}".format(y_k[i] * tmp1, (2-y_k[i]) * tmp1)) + feed_back = "integrand nan" + break + r += r_tmp + I.append(tmp1*r) + + return np.asarray(I), feed_back + + + def get_fourier_integral_simps_weighted_values(yl): N = len(yl) if N % 2 == 1: # odd N @@ -338,7 +397,7 @@ def get_dt_for_accurate_interpolation(t_max, tol, ft_ref, diff_method=_absDiff): N*=2 -def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_only, N_max = 2**20, diff_method=_absDiff): +def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, ft_ref, opt_b_only, diff_method=_absDiff): log.info("get_dt_for_accurate_interpolation, please wait ...") try: @@ -347,8 +406,6 @@ def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_ except RuntimeError: c = t_max - - c = min(c, t_max) dt_tol = get_dt_for_accurate_interpolation(t_max=c, tol=intpl_tol, @@ -369,27 +426,25 @@ def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_ dt = 2*np.pi/dx/N - if dt <= dt_tol: - log.debug("dt criterion fulfilled") - return a, b, N, dx, dt - else: + if dt > dt_tol: log.debug("down scale dx and dt to match new power of 2 N") - N_min = 2*np.pi/dx/dt_tol - N = 2**int(np.ceil(np.log2(N_min))) - scale = np.sqrt(N_min/N) - dx_new = scale*dx - b_minus_a = dx_new*N - dt_new = 2*np.pi/dx_new/N - assert dt_new < dt_tol + N_min = 2*np.pi/dx/dt_tol + N = 2**int(np.ceil(np.log2(N_min))) + scale = np.sqrt(N_min/N) + dx_new = scale*dx + b_minus_a = dx_new*N + dt_new = 2*np.pi/dx_new/N + assert dt_new < dt_tol - print(a, b) - if opt_b_only: - b = a + b_minus_a + if opt_b_only: + b = a + b_minus_a + else: + delta = b_minus_a - (b-a) + b += delta/2 + a -= delta/2 else: - delta = b_minus_a - (b-a) - b += delta/2 - a -= delta/2 + dt_new = dt if dt_new*(N-1) < t_max: log.info("increase N to match dt_new*(N-1) < t_max") @@ -397,7 +452,4 @@ def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_ N = 2 ** int(np.ceil(np.log2(N_tmp))) dx_new = 2*np.pi/N/dt_new - print("dx", dx_new) - print("(b-a)/N", (b-a)/N) - return a, b, N, dx_new, dt_new diff --git a/stocproc/stocproc.py b/stocproc/stocproc.py index cab8ff9..6588812 100644 --- a/stocproc/stocproc.py +++ b/stocproc/stocproc.py @@ -357,51 +357,29 @@ class StocProc_FFT(_absStocProc): if not negative_frequencies: log.info("non neg freq only") - # assume the spectral_density is 0 for w<0 - # and decays fast for large w - b = method_fft.find_integral_boundary(integrand = spectral_density, - tol = intgr_tol**2, - ref_val = 1, - max_val = 1e6, - x0 = 0.777) - log.info("upper int bound b {:.3e}".format(b)) a, b, N, dx, dt = method_fft.calc_ab_N_dx_dt(integrand = spectral_density, intgr_tol = intgr_tol, intpl_tol = intpl_tol, t_max = t_max, - a = 0, - b = b, ft_ref = lambda tau:bcf_ref(tau)*np.pi, - opt_b_only= True, - N_max = 2**24) - log.info("required tol results in N {}".format(N)) + opt_b_only= True) else: log.info("use neg freq") - # assume the spectral_density is non zero also for w<0 - # but decays fast for large |w| - b = method_fft.find_integral_boundary(integrand = spectral_density, - tol = intgr_tol, - ref_val = 1, - max_val = 1e6, - x0 = 0.777) - a = method_fft.find_integral_boundary(integrand = spectral_density, - tol = intgr_tol, - ref_val = -1, - max_val = 1e6, - x0 = -0.777) a, b, N, dx, dt = method_fft.calc_ab_N_dx_dt(integrand = spectral_density, intgr_tol = intgr_tol, intpl_tol = intpl_tol, t_max = t_max, - a = a, - b = b, ft_ref = lambda tau:bcf_ref(tau)*np.pi, - opt_b_only= False, - N_max = 2**24) - log.info("required tol result in N {}".format(N)) + opt_b_only= False) 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 @@ -447,4 +425,127 @@ class StocProc_FFT(_absStocProc): r"""The number of independent random variables Y is given by the number of discrete times :math:`t_l < t_\mathrm{max}` from the Fourier Transform """ - return len(self.yl) \ No newline at end of file + return len(self.yl) + + +class StocProc_TanhSinh(_absStocProc): + r"""Simulate Stochastic Process using TanhSinh integration for the Fourier Integral + """ + + def __init__(self, spectral_density, t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2, + seed=None, negative_frequencies=False, scale=1): + self.key = bcf_ref, t_max, intgr_tol, intpl_tol + + if not negative_frequencies: + log.info("non neg freq only") + log.info("get_dt_for_accurate_interpolation, please wait ...") + try: + ft_ref = lambda tau: bcf_ref(tau) * np.pi + c = method_fft.find_integral_boundary(lambda tau: np.abs(ft_ref(tau)) / np.abs(ft_ref(0)), + intgr_tol, 1, 1e6, 0.777) + except RuntimeError: + c = t_max + + c = min(c, t_max) + dt_tol = method_fft.get_dt_for_accurate_interpolation(t_max=c, + tol=intpl_tol, + ft_ref=ft_ref) + log.info("requires dt < {:.3e}".format(dt_tol)) + else: + raise NotImplementedError + + N = int(np.ceil(t_max/dt_tol)) + log.info("yields N = {}".format(N)) + + wmax = method_fft.find_integral_boundary(spectral_density, tol=intgr_tol/4, ref_val=1, max_val=1e6, x0=0.777) + + h = 0.1 + k = 15 + kstep = 5 + conv_fac = 2 + old_d = None + log.info("find h and kmax for TanhSinh integration ...") + while True: + tau = np.asarray([t_max]) + num_FT, fb = method_fft.fourier_integral_TanhSinh_with_feedback(integrand=lambda w: spectral_density(w) / np.pi, + w_max=wmax, + tau=tau, + h=h, + kmax=k) + d = np.abs(num_FT - bcf_ref(tau)) / np.abs(bcf_ref(0)) + 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)) + h /= 2 + k = 15 + kstep *= 2 + + if old_d is None: + old_d = d[0] + else: + if old_d < conv_fac * d[0]: + wmax *= 1.5 + log.info("convergence factor of {} not met -> inc wmax to {}".format(conv_fac, wmax)) + h = 0.1 + k = 15 + kstep = 5 + old_d = None + else: + old_d = d[0] + + if d < intgr_tol: + log.info("intgration tolerance met with h {} and kmax {}".format(h, k)) + break + + tau = np.linspace(0, (N-1)*dt_tol, N) + num_FT = method_fft.fourier_integral_TanhSinh( + integrand=lambda w: spectral_density(w) / np.pi, + w_max=wmax, + tau=tau, + h=h, + kmax=k) + d = np.max(np.abs(num_FT - bcf_ref(tau)) / np.abs(bcf_ref(0))) + assert d < intgr_tol + + 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)] + + tmp1 = wmax / 2 + #sd = np.hstack((spectral_density(yk * tmp1, )) + + + + + self.yl = wk*spectral_density() + + + + + + + + # 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