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

unmerged changes from some late night work

Browse source
This commit is contained in:
Richard Hartmann 2016-10-27 22:51:48 +02:00
parent cfbf5fffc8
commit b7c3908f06
2 changed files with 129 additions and 70 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

69
stocproc/method_fft.py
View file

@ -1,8 +1,10 @@
from scipy.optimize import brentq from scipy.optimize import brentq, bisect
from numpy.fft import rfft as np_rfft from numpy.fft import rfft as np_rfft
import numpy as np 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): def find_integral_boundary(integrand, tol, ref_val, max_val, x0):
""" """
@ -23,23 +25,45 @@ def find_integral_boundary(integrand, tol, ref_val, max_val, x0):
if integrand(ref_val) <= tol: if integrand(ref_val) <= tol:
raise ValueError("the integrand at ref_val needs to be greater that tol") raise ValueError("the integrand at ref_val needs to be greater that tol")
# find the left boundary called a log.debug("ref_value for search: {} tol: {}".format(ref_val, tol))
if integrand(x0+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 x = 2*x0
while integrand(x+ref_val) > tol: I = integrand(x + ref_val)
if abs(x-ref_val) > max_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") 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) 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 x = x0/2
while integrand(x+ref_val) < tol: I = integrand(x + ref_val)
if (1/abs(x-ref_val)) > max_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") raise RuntimeError("1/|x-ref_val| > max_val was reached")
x /= 2 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: else:
a = x0 a = x0
log.debug("I(ref_val) = tol -> no search necessary")
log.debug("return a={:.5g}".format(a))
return a return a
def find_integral_boundary_auto(integrand, tol, ref_val=0, max_val=1e6, def find_integral_boundary_auto(integrand, tol, ref_val=0, max_val=1e6,
@ -51,44 +75,49 @@ def find_integral_boundary_auto(integrand, tol, ref_val=0, max_val=1e6,
max_val_left = max_val if max_val_left is None else max_val_left 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 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) 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) b = find_integral_boundary(integrand, tol, ref_val=ref_val_right, max_val=max_val_right, x0=+1)
return a,b return a,b
def fourier_integral(integrand, a, b, N): def fourier_integral(integrand, a, b, N):
""" """
approximates int_a^b dx integrand(x) by the riemann sum with N terms 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_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)) yl = integrand(np.linspace(a+delta_x/2, b+delta_x/2, N, endpoint=False))
fft_vals = np_rfft(yl) fft_vals = np_rfft(yl)
tau = np.arange(len(fft_vals))*delta_k 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 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 chooses N such that the approximated Fourier integral
meets the exact solution within a given tolerance of the meets the exact solution within a given tolerance of the
relative deviation for a given interval of interest 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 i = 10
while True: while True:
N = 2**i N = 2**i
tau, ft_tau = fourier_integral(integrand, a, b, N) tau, ft_tau = fourier_integral(integrand, a, b, N)
idx = np.where(tau <= t_max) idx = np.where(tau <= t_max)
ft_ref_tau = ft_ref(tau[idx]) 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: if rd < tol:
log.debug("reached rd ({:.3e}) < tol ({:.3e}), return N={}".format(rd, tol, N))
return N return N
if N > N_max:
raise RuntimeError("maximum number of points for Fourier Transform reached")
i += 2 i += 1

106
tests/test_method_fft.py
View file

@ -12,6 +12,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)
@ -78,6 +83,25 @@ def fourier_integral_trapz(integrand, a, b, N):
return tau, delta_x*np.exp(-1j*tau*a)*fft_vals 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): def fourier_integral_simple_test(integrand, a, b, N):
delta_x = (b-a)/N delta_x = (b-a)/N
delta_k = 2*np.pi/(b-a) delta_k = 2*np.pi/(b-a)
@ -119,6 +143,23 @@ def fourier_integral_trapz_simple_test(integrand, a, b, N):
return k, res 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(): def test_fourier_integral():
intg = lambda x: x**2 intg = lambda x: x**2
a = -1.23 a = -1.23
@ -149,47 +190,28 @@ def test_fourier_integral():
assert rd < 4e-6 assert rd < 4e-6
N = 512 N = 1024
tau, ft_n = sp.method_fft.fourier_integral(intg, a, b, N) for fac in [1,2,4,8]:
ft_ref_n = ft_ref(tau) 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_simple = fourier_integral_simple_test(intg, a, b, N) k, ft_n = fourier_integral_trapz(intg, a, b, fac*N)
assert np.max(np.abs(ft_simple-ft_n)) < 1e-11 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))
# plt.plot(np.abs(ft_simple-ft_n)) k, ft_n = fourier_integral_simps(intg, a, b, fac*N-1)
# plt.yscale('log') ft_ref_n = ft_ref(k)
# plt.show() # 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()
N = 512 plt.yscale('log')
tau, ft_n = fourier_integral_trapz(intg, a, b, N) plt.legend()
ft_ref_n = ft_ref(tau) plt.show()
sys.exit()
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 s = 0.1
wc = 4 wc = 4
@ -225,11 +247,19 @@ def test_fourier_integral():
# plt.yscale('log') # plt.yscale('log')
# plt.show() # 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__": if __name__ == "__main__":
# test_find_integral_boundary() # test_find_integral_boundary()
test_fourier_integral() test_fourier_integral()
# test_get_N_for_accurate_fourier_integral()