import os import pathlib import sys p = pathlib.PosixPath(os.path.abspath(__file__)) sys.path.insert(0, str(p.parent.parent)) import math import logging import numpy as np import scipy.integrate as sp_int from scipy.special import gamma as gamma_func import stocproc as sp from stocproc import method_fft try: import matplotlib.pyplot as plt except ImportError: print("matplotlib not found -> any plotting will crash") def test_find_integral_boundary(): def f(x): return np.exp(-(x)**2) tol = 1e-10 b = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=+1, max_val=1e6) a = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=-1, max_val=1e6) assert a != b assert abs(f(a)-tol) < 1e-14 assert abs(f(b)-tol) < 1e-14 b = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=b+5, max_val=1e6) a = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=a-5, max_val=1e6) assert a != b assert abs(f(a)-tol) < 1e-14 assert abs(f(b)-tol) < 1e-14 def f2(x): return np.exp(-(x)**2)*x**2 tol = 1e-10 b = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=1, x0=+1, max_val=1e6) a = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=-1, x0=-1, max_val=1e6) assert a != b assert abs(f2(a) -tol) < 1e-14 assert abs(f2(b)-tol) < 1e-14 b = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=1, x0=b+5, max_val=1e6) a = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=-1, x0=a-5, max_val=1e6) assert a != b assert abs(f2(a)-tol) < 1e-14, "diff {}".format(abs(f2(a)/f2(1)-tol)) assert abs(f2(b)-tol) < 1e-14, "diff {}".format(abs(f2(b)/f2(-1)-tol)) def f3(x): return np.exp(-(x-5)**2)*x**2 tol = 1e-10 b = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=+1, max_val=1e6) a = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=-1, max_val=1e6) assert a != b assert abs(f3(a)-tol) < 1e-14 assert abs(f3(b)-tol) < 1e-14 b = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=b+5, max_val=1e6) a = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=a-5, max_val=1e6) assert a != b assert abs(f3(a)-tol) < 1e-14, "diff {}".format(abs(f3(a)-tol)) assert abs(f3(b)-tol) < 1e-14, "diff {}".format(abs(f3(b)-tol)) ################################## ## the case where f(xref) < tol ## ################################## def f(x): return np.exp(-x ** 2) tol = 1e-3 b = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=10, x0=+1, max_val=1e6) assert abs(f(b) - tol) < 1e-14 a = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=-10, x0=-1., max_val=1e6) assert abs(f(a) - tol) < 1e-14 def fourier_integral_trapz(integrand, a, b, N): """ approximates int_a^b dx integrand(x) by the riemann sum with N terms this function is here and not in method_fft because it has almost no advantage over the modpoint method. so only for testing purposes. """ yl = integrand(np.linspace(a, b, N+1, endpoint=True)) yl[0] = yl[0]/2 yl[-1] = yl[-1]/2 delta_x = (b-a)/N delta_k = 2*np.pi*N/(b-a)/(N+1) fft_vals = np.fft.rfft(yl) tau = np.arange(len(fft_vals))*delta_k return tau, delta_x*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) x = np.linspace(a, b, N, endpoint = False) + delta_x/2 k = np.arange(N//2+1)*delta_k k_np = 2*np.pi*np.fft.rfftfreq(N, delta_x) kdif = np.max(np.abs(k_np[1:] - k[1:])/k[1:]) if kdif > 1e-15: print("WARNING |rfftfreq - k| = {}".format(kdif)) yl = integrand(x) res = np.empty(shape=(N//2+1,), dtype=np.complex128) for i in range(N//2+1): tmp = yl*np.exp(-1j*x*k[i]) res[i] = delta_x*(math.fsum(tmp.real) + 1j*math.fsum(tmp.imag)) return k, res def fourier_integral_trapz_simple_test(integrand, a, b, N): delta_x = (b-a)/N delta_k = 2*np.pi*N/(b-a)/(N+1) x = np.linspace(a, b, N+1, endpoint = True) k = np.arange((N+1)//2+1)*delta_k yl = integrand(x) yl[0] = yl[0]/2 yl[-1] = yl[-1]/2 res = np.empty(shape=((N+1)//2+1,), dtype=np.complex128) for i in range((N+1)//2+1): tmp = yl*np.exp(-1j*x*k[i]) res[i] = delta_x*(math.fsum(tmp.real) + 1j*math.fsum(tmp.imag)) return k, res def test_fourier_integral_finite_boundary(): intg = lambda x: x**2 a = -1.23 b = 4.87 ################################ ## check with analytic result ## ################################ ft_ref = lambda k: (np.exp(-1j*a*k)*(2j - a*k*(2 + 1j*a*k)) + np.exp(-1j*b*k)*(-2j + b*k*(2 + 1j*b*k)))/k**3 N = 2**18 N_test = 100 tau, ft_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) ft_ref_n = ft_ref(tau) tau = tau[1:N_test] ft_n = ft_n[1:N_test] ft_ref_n = ft_ref_n[1:N_test] rd = np.max(np.abs(ft_n - ft_ref_n)/np.abs(ft_ref_n)) assert rd < 4e-6, "rd = {}".format(rd) N = 2**18 N_test = 100 tau, ft_n = fourier_integral_trapz(intg, a, b, N) ft_ref_n = ft_ref(tau) tau = tau[1:N_test] ft_n = ft_n[1:N_test] ft_ref_n = ft_ref_n[1:N_test] rd = np.max(np.abs(ft_n - ft_ref_n)/np.abs(ft_ref_n)) assert rd < 4e-6, "rd = {}".format(rd) ###################################################### ## check against numeric fourier integral (non FFT) ## ###################################################### N = 512 tau, ft_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) k, ft_simple = fourier_integral_simple_test(intg, a, b, N) assert np.max(np.abs(ft_simple-ft_n)) < 1e-11 N = 512 tau, ft_n = fourier_integral_trapz(intg, a, b, N) k, ft_simple = fourier_integral_trapz_simple_test(intg, a, b, N) assert np.max(np.abs(ft_simple-ft_n)) < 1e-11 ################################# ## check midp against simpson ## ################################# N = 1024 tau, ft_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) ft_ref_n = ft_ref(tau) rd = np.abs(ft_ref_n-ft_n) / np.abs(ft_ref_n) idx = np.where(np.logical_and(tau < 75, np.isfinite(rd))) rd = rd[idx] mrd_midp = np.max(rd) assert mrd_midp < 9e-3, "mrd_midp = {}".format(mrd_midp) N = 513 tau, ft_n = sp.method_fft.fourier_integral_simps(intg, a, b, N) ft_ref_n = ft_ref(tau) rd = np.abs(ft_ref_n-ft_n) / np.abs(ft_ref_n) idx = np.where(np.logical_and(tau < 75, np.isfinite(rd))) rd = rd[idx] mrd_simps = np.max(rd) assert mrd_simps < 4e-3, "mrd_simps = {}".format(mrd_midp) assert mrd_simps < mrd_midp, "mrd_simps ({:.3e}) >= mrd_trapz ({:.3e})".format(mrd_simps, mrd_trapz) 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[idx_pos]**s*np.exp(-w_flat[idx_pos]/wc) return res def test_fourier_integral_infinite_boundary(plot=False): s = 0.5 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) errs = [9e-5, 2e-5, 1.3e-6] for i, N in enumerate([2**16, 2**18, 2**20]): tau, bcf_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N=N) bcf_ref_n = bcf_ref(tau) tau_max = 5 idx = np.where(tau <= tau_max) tau = tau[idx] bcf_n = bcf_n[idx] bcf_ref_n = bcf_ref_n[idx] rd_mp = np.abs(bcf_ref_n-bcf_n)/np.abs(bcf_ref_n) if plot: p, = plt.plot(tau, rd_mp, label="trapz N {}".format(N)) tau, bcf_n = sp.method_fft.fourier_integral_simps(intg, a, b=b, N=N-1) bcf_ref_n = bcf_ref(tau) idx = np.where(tau <= tau_max) tau = tau[idx] bcf_n = bcf_n[idx] bcf_ref_n = bcf_ref_n[idx] rd_sm = np.abs(bcf_ref_n-bcf_n)/np.abs(bcf_ref_n) if plot: plt.plot(tau, rd_sm, label="simps N {}".format(N), color=p.get_color(), ls='--') t_ = 3 x_simps, dx = np.linspace(a,b,N-1, endpoint=True, retstep=True) I = sp_int.simps(intg(x_simps)*np.exp(-1j*x_simps*t_), dx=dx) err = np.abs(I-bcf_ref(t_))/np.abs(bcf_ref(t_)) assert np.max(rd_mp) < errs[i], "np.max(rd_mp) = {} >= {}".format(np.max(rd_mp), errs[i]) assert np.max(rd_sm) < errs[i], "np.max(rd_sm) = {} >= {}".format(np.max(rd_sm), errs[i]) if plot: plt.plot(t_, err, marker='o', color='g') if plot: plt.legend(loc='lower right') plt.grid() plt.yscale('log') plt.show() def test_get_N_a_b_for_accurate_fourier_integral(): _WC_ = 2 intg = lambda w: 1 / (1 + (w - _WC_) ** 2) / np.pi bcf_ref = lambda t: np.exp(- np.abs(t) - 1j * _WC_ * t) a, b = sp.method_fft.find_integral_boundary_auto(integrand=intg, tol=1e-2, ref_val=_WC_) N, a, b = sp.method_fft.get_N_a_b_for_accurate_fourier_integral(intg, t_max=50, tol=1e-2, ft_ref=bcf_ref, opt_b_only=False) print(N, a, b) def test_get_N_a_b_for_accurate_fourier_integral_b_only(): s = 0.5 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-2, ref_val=1) a = 0 N, a, b = sp.method_fft.get_N_a_b_for_accurate_fourier_integral(intg, t_max=15, tol=1e-5, ft_ref=bcf_ref, opt_b_only=True) print(N,a,b) def test_get_dt_for_accurate_interpolation(): s = 0.5 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)) dt = sp.method_fft.get_dt_for_accurate_interpolation(t_max=40, tol=1e-4, ft_ref=bcf_ref) print(dt) def test_sclicing(): yl = np.ones(10, dtype=int) yl = sp.method_fft.get_fourier_integral_simps_weighted_values(yl) assert yl[0] == 2/6 assert yl[1] == 8/6 assert yl[2] == 4/6 assert yl[3] == 8/6 assert yl[4] == 4/6 assert yl[5] == 8/6 assert yl[6] == 4/6 assert yl[7] == 8/6 assert yl[8] == 5/6 assert yl[9] == 3/6 def test_calc_abN(): def testing(intg, bcf_ref, tol, tmax): diff_method = method_fft._absDiff a, b = sp.method_fft.find_integral_boundary_auto(integrand=intg, tol=1e-2, ref_val=1) a, b, N, dx, dt = sp.method_fft.calc_ab_N_dx_dt(integrand=intg, intgr_tol=tol, intpl_tol=tol, t_max=tmax, a=0, b=b, ft_ref=bcf_ref, opt_b_only=True, diff_method=diff_method) tau, ft_tau = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) idx = np.where(tau <= tmax) ft_ref_tau = bcf_ref(tau[idx]) rd = diff_method(ft_tau[idx], ft_ref_tau) tau_fine = np.linspace(0, tmax, 1500) ft_ref_n = bcf_ref(tau_fine) ft_intp = sp.tools.ComplexInterpolatedUnivariateSpline(x=tau[idx], y=ft_tau[idx], k=3, noWarning=True) ft_intp_n = ft_intp(tau_fine) d = diff_method(ft_intp_n, ft_ref_n) assert rd < tol assert d < tol assert (np.abs(dx * dt * N - np.pi * 2)) < 1e-15 s = 0.5 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)) tol = 1e-3 tmax=40 testing(intg, bcf_ref, tol, tmax) s = 0.5 wc = 40 intg = lambda x: osd(x, s, wc) bcf_ref = lambda t: gamma_func(s + 1) * wc ** (s + 1) * (1 + 1j * wc * t) ** (-(s + 1)) tol = 1e-3 tmax = 40 testing(intg, bcf_ref, tol, tmax) 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] def sd(w): s = -0.5 wc = 5 return w**s*np.exp(-w/wc) def bcf(tau): s = -0.5 wc = 5 return (wc/(1+1j*wc*tau))**(s+1)*gamma_func(s+1)/np.pi def test_SP_TanhSinh_dump(): tmax = 25 _sp = sp.StocProc_TanhSinh(spectral_density=sd, t_max=tmax, bcf_ref=bcf, intgr_tol=1e-2, intpl_tol=1e-2, seed=0) t = np.linspace(0, tmax, 500) _sp.new_process() zt = _sp(t) import pickle sp_dump = pickle.dumps(_sp) del _sp _sp2 = pickle.loads(sp_dump) _sp2.new_process() zt2 = _sp2(t) assert np.all(zt == zt2) N = 1000 t = np.linspace(0, tmax, 500) idx = 200 c = 0 for i in range(N): _sp2.new_process() zt = _sp2(t) c += zt*np.conj(zt[idx]) c /= N d = np.max(np.abs(c - bcf(t-t[idx]))) print(d) assert d < 0.09 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() # test_SP_TanhSinh() test_SP_TanhSinh_dump()