mirror of
https://github.com/vale981/stocproc
synced 2025-03-04 17:21:42 -05:00
work in progress on TanhSinh SP for singular SD
This commit is contained in:
Richard Hartmann
parent
9a061f5426
commit
65f4fb4eee
2 changed files with 207 additions and 54 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
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
|
||||
Loading…
Add table
Reference in a new issue