From b7c3908f06c5526c75ea8184868d8f2387baa5f4 Mon Sep 17 00:00:00 2001 From: Richard Hartmann Date: Thu, 27 Oct 2016 22:51:48 +0200 Subject: [PATCH] unmerged changes from some late night work --- stocproc/method_fft.py | 77 ++++++++++++++++-------- tests/test_method_fft.py | 122 ++++++++++++++++++++++++--------------- 2 files changed, 129 insertions(+), 70 deletions(-) diff --git a/stocproc/method_fft.py b/stocproc/method_fft.py index 35f24aa..9cba80f 100644 --- a/stocproc/method_fft.py +++ b/stocproc/method_fft.py @@ -1,8 +1,10 @@ -from scipy.optimize import brentq +from scipy.optimize import brentq, bisect from numpy.fft import rfft as np_rfft import numpy as np -from math import fsum +import logging + +log = logging.getLogger(__name__) def find_integral_boundary(integrand, tol, ref_val, max_val, x0): """ @@ -22,24 +24,46 @@ def find_integral_boundary(integrand, tol, ref_val, max_val, x0): assert x0 != 0 if integrand(ref_val) <= tol: raise ValueError("the integrand at ref_val needs to be greater that tol") - - # find the left boundary called a - if integrand(x0+ref_val) > tol: + + log.debug("ref_value for search: {} tol: {}".format(ref_val, tol)) + + I = integrand(x0+ref_val) + + if I > tol: + log.debug("x={:.3e} I(x+ref_val) = {:.3e} > tol -> veer x away from ref_value".format(x0, I)) x = 2*x0 - while integrand(x+ref_val) > tol: - if abs(x-ref_val) > max_val: + I = integrand(x + ref_val) + while I > tol: + log.debug("x={:.3e} I(x+ref_val) = {:.3e}".format(x, I)) + if abs(x) > max_val: raise RuntimeError("|x-ref_val| > max_val was reached") - x *= 2 + x *= 2 + I = integrand(x + ref_val) + + log.debug("x={:.3e} I(x+ref_val) = {:.3e} < tol".format(x, I)) a = brentq(lambda x: integrand(x)-tol, x+ref_val, x0+ref_val) - elif integrand(x0+ref_val) < tol: + log.debug("found I(a={:.3e}) = {:.3e} = tol".format(a, integrand(a))) + + elif I < tol: + log.debug("x={:.3e} I(x+ref_val) = {:.3e} < tol -> approach x towards ref_value".format(x0, I)) x = x0/2 - while integrand(x+ref_val) < tol: - if (1/abs(x-ref_val)) > max_val: + I = integrand(x + ref_val) + while I < tol: + log.debug("x={:.3e} I(x+ref_val) = {:.3e}".format(x, I)) + if (1/abs(x)) > max_val: raise RuntimeError("1/|x-ref_val| > max_val was reached") x /= 2 - a = brentq(lambda x: integrand(x)-tol, x+ref_val, x0+ref_val) + I = integrand(x+ref_val) + + log.debug("x={:.3e} I(x+ref_val) = {:.3e} > tol".format(x, I)) + log.debug("search for root in interval [{:.3e},{:.3e}]".format(x0+ref_val, x+ref_val)) + a = brentq(lambda x_: integrand(x_)-tol, x+ref_val, x0+ref_val) + log.debug("found I(a={:.3e}) = {:.3e} = tol".format(a, integrand(a))) else: a = x0 + log.debug("I(ref_val) = tol -> no search necessary") + + log.debug("return a={:.5g}".format(a)) return a def find_integral_boundary_auto(integrand, tol, ref_val=0, max_val=1e6, @@ -50,45 +74,50 @@ def find_integral_boundary_auto(integrand, tol, ref_val=0, max_val=1e6, ref_val_right = ref_val if ref_val_right is None else ref_val_right max_val_left = max_val if max_val_left is None else max_val_left max_val_right = max_val if max_val_right is None else max_val_right - + + log.debug("trigger left search") a = find_integral_boundary(integrand, tol, ref_val=ref_val_left, max_val=max_val_left, x0=-1) + log.debug("trigger right search") b = find_integral_boundary(integrand, tol, ref_val=ref_val_right, max_val=max_val_right, x0=+1) return a,b def fourier_integral(integrand, a, b, N): """ approximates int_a^b dx integrand(x) by the riemann sum with N terms - + and the most simplest uniform midpoint weights """ + log.debug("integrate over [{:.3e},{:.3e}] using {} points".format(a,b,N)) delta_x = (b-a)/N - delta_k = 2*np.pi/(b-a) + delta_k = 2*np.pi/(b-a) yl = integrand(np.linspace(a+delta_x/2, b+delta_x/2, N, endpoint=False)) fft_vals = np_rfft(yl) tau = np.arange(len(fft_vals))*delta_k + 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 get_N_for_accurate_fourier_integral(integrand, a, b, t_0, t_max, tol, ft_ref, N_max = 2**15): +def get_N_for_accurate_fourier_integral(integrand, a, b, t_max, ft_ref, tol=1e-3, N_max = 2**20): """ chooses N such that the approximated Fourier integral meets the exact solution within a given tolerance of the relative deviation for a given interval of interest """ + + log.debug("FFT accuracy for k in [0, {:.3e}] (tol={:.3e})".format(t_max, tol)) + i = 10 while True: N = 2**i tau, ft_tau = fourier_integral(integrand, a, b, N) idx = np.where(tau <= t_max) ft_ref_tau = ft_ref(tau[idx]) - rd = np.abs(ft_tau[idx] - ft_ref_tau) / np.abs(ft_ref_tau) + rd = np.max(np.abs(ft_tau[idx] - ft_ref_tau) / np.abs(ft_ref_tau)) + log.debug("N:{} found rel dif of {:.3e}".format(N, rd)) if rd < tol: + log.debug("reached rd ({:.3e}) < tol ({:.3e}), return N={}".format(rd, tol, N)) return N + if N > N_max: + raise RuntimeError("maximum number of points for Fourier Transform reached") - i += 2 + i += 1 - - - - - - \ No newline at end of file diff --git a/tests/test_method_fft.py b/tests/test_method_fft.py index b9bf221..8558a1a 100644 --- a/tests/test_method_fft.py +++ b/tests/test_method_fft.py @@ -12,6 +12,11 @@ sys.path.insert(0, str(p.parent.parent)) import stocproc as sp +import logging +logging.basicConfig(level=logging.DEBUG) + + + def test_find_integral_boundary(): def f(x): return np.exp(-(x)**2) @@ -78,6 +83,25 @@ def fourier_integral_trapz(integrand, a, b, N): return tau, delta_x*np.exp(-1j*tau*a)*fft_vals + +def fourier_integral_simps(integrand, a, b, N): + """ + approximates int_a^b dx integrand(x) by the riemann sum with N terms + + """ + assert (N % 2) == 1 + yl = integrand(np.linspace(a, b, N, endpoint=True)) + yl[1::2] *= 4 + yl[2:-2:2] *= 2 + + delta_x = (b - a) / (N-1) + delta_k = 2 * np.pi / (N*delta_x) + + fft_vals = np.fft.rfft(yl) + tau = np.arange(len(fft_vals)) * delta_k + + return tau, delta_x/3 * np.exp(-1j * tau * a) * fft_vals + def fourier_integral_simple_test(integrand, a, b, N): delta_x = (b-a)/N delta_k = 2*np.pi/(b-a) @@ -118,6 +142,23 @@ def fourier_integral_trapz_simple_test(integrand, a, b, N): res[i] = delta_x*(math.fsum(tmp.real) + 1j*math.fsum(tmp.imag)) return k, res + + +def osd(w, s, wc): + if not isinstance(w, np.ndarray): + if w < 0: + return 0 + else: + return w ** s * np.exp(-w / wc) + else: + res = np.zeros(shape=w.shape) + + w_flat = w.flatten() + idx_pos = np.where(w_flat > 0) + fv_res = res.flat + fv_res[idx_pos] = w_flat ** s * np.exp(-w_flat / wc) + + return res def test_fourier_integral(): intg = lambda x: x**2 @@ -148,49 +189,30 @@ def test_fourier_integral(): # print(rd) assert rd < 4e-6 + + N = 1024 + for fac in [1,2,4,8]: + k, ft_n = sp.method_fft.fourier_integral(intg, a, b, fac*N) + ft_ref_n = ft_ref(k) + # assert np.max(np.abs(ft_n-ft_ref_n)) < 1e-11 + plt.plot(k, np.abs(ft_n-ft_ref_n)/np.abs(ft_ref_n), label='simple f:{}'.format(fac)) + + k, ft_n = fourier_integral_trapz(intg, a, b, fac*N) + ft_ref_n = ft_ref(k) + # assert np.max(np.abs(ft_n-ft_ref_n)) < 1e-11 + plt.plot(k, np.abs(ft_n-ft_ref_n)/np.abs(ft_ref_n), label='trapz f:{}'.format(fac)) + + k, ft_n = fourier_integral_simps(intg, a, b, fac*N-1) + ft_ref_n = ft_ref(k) + # assert np.max(np.abs(ft_n-ft_ref_n)) < 1e-11 + plt.plot(k, np.abs(ft_n-ft_ref_n)/np.abs(ft_ref_n), label='simps f:{}'.format(fac)) + + plt.grid() + plt.yscale('log') + plt.legend() + plt.show() + sys.exit() - N = 512 - tau, ft_n = sp.method_fft.fourier_integral(intg, a, b, N) - ft_ref_n = ft_ref(tau) - - k, ft_simple = fourier_integral_simple_test(intg, a, b, N) - assert np.max(np.abs(ft_simple-ft_n)) < 1e-11 - -# plt.plot(np.abs(ft_simple-ft_n)) -# plt.yscale('log') -# plt.show() - - - N = 512 - tau, ft_n = fourier_integral_trapz(intg, a, b, N) - ft_ref_n = ft_ref(tau) - - k, ft_simple = fourier_integral_trapz_simple_test(intg, a, b, N) - assert np.max(np.abs(ft_simple-ft_n)) < 1e-11 - -# plt.plot(np.abs(ft_simple-ft_n)) -# plt.yscale('log') -# plt.show() -# sys.exit() - - - - def osd(w, s, wc): - if not isinstance(w, np.ndarray): - if w < 0: - return 0 - else: - return w**s*np.exp(-w/wc) - else: - res = np.zeros(shape=w.shape) - - w_flat = w.flatten() - idx_pos = np.where(w_flat > 0) - fv_res = res.flat - fv_res[idx_pos] = w_flat**s*np.exp(-w_flat/wc) - - return res - s = 0.1 wc = 4 @@ -225,11 +247,19 @@ def test_fourier_integral(): # plt.yscale('log') # plt.show() - - +def test_get_N_for_accurate_fourier_integral(): + s = 0.1 + wc = 4 + + intg = lambda x: osd(x, s, wc) + 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) + + N = sp.method_fft.get_N_for_accurate_fourier_integral(intg, a, b, t_max = 40, ft_ref=bcf_ref, tol = 1e-3, N_max = 2**20) if __name__ == "__main__": -# test_find_integral_boundary() - test_fourier_integral() \ No newline at end of file + # test_find_integral_boundary() + test_fourier_integral() + # test_get_N_for_accurate_fourier_integral()