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

all tests pass

Browse source
This commit is contained in:
Richard Hartmann 2016-11-15 14:36:33 +01:00
parent d58c0fde1a
commit b3b7b1064d
4 changed files with 118 additions and 81 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

89
stocproc/method_fft.py
View file

@ -43,14 +43,43 @@ def find_integral_boundary(integrand, tol, ref_val, max_val, x0):
"""
_max_num_iteration = 100
_i = 0
x0 = float(x0)
I_ref = integrand(ref_val)
if I_ref < tol:
pass
# return find_integral_boundary(integrand = lambda x: 1/integrand(x),
# tol = 1/tol,
# ref_val=ref_val,
# max_val=max_val,
# x0=-x0)
x_old = ref_val
while True:
_i += 1
if _i > _max_num_iteration:
raise RuntimeError("max number of iteration reached")
x = ref_val - x0
try:
I_x = integrand(x)
except Exception as e:
raise RuntimeError("evaluation of integrand failed due to {}".format(e))
if I_x > tol:
break
x0 *= 2
x_old = x
a = brentq(lambda x: integrand(x) - tol, x_old, x)
return a
elif I_ref > tol:
x_old = ref_val
while True:
_i += 1
if _i > _max_num_iteration:
raise RuntimeError("max number of iteration reached")
x = ref_val + x0
I_x = integrand(x)
try:
I_x = integrand(x)
except Exception as e:
raise RuntimeError("evaluation of integrand failed due to {}".format(e))
if I_x < tol:
break
x0 *= 2
@ -148,28 +177,33 @@ def _f_opt(x, integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, b_on
else:
a_ = find_integral_boundary(integrand, tol=tol, ref_val=a, max_val=1e6, x0=-1)
b_ = find_integral_boundary(integrand, tol=tol, ref_val=b, max_val=1e6, x0=1)
except ValueError:
a_ = -1
b_ = 1
except RuntimeError:
# in case 'find_integral_boundary' failes
d = 300
_f_opt_cache[key] = d, None, None
return d
tau, ft_tau = fourier_integral_midpoint(integrand, a_, b_, N)
idx = np.where(tau <= t_max)
ft_ref_tau = ft_ref(tau[idx])
d = diff_method(ft_ref_tau, ft_tau[idx])
_f_opt_cache[key] = d, a_, b_
#log.debug("f_opt tol {} -> d {}".format(10**x, d))
return np.log10(d)
def _lower_contrs(x, integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, b_only):
_f_opt(x, integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, b_only)
tol = 10**x
d, a_, b_ = _f_opt_cache[float(x[0])]
if (a_ is None) or (b_ is None):
return -1
v = N * np.pi / (b_ - a_) - t_max
log.debug("lower constr value {} for x {} (tol {})".format(v, x, tol))
#log.debug("lower constr value {} for x {} (tol {})".format(v, 10**x, tol))
return v
def _upper_contrs(x):
log.debug("upper constr value {}".format(-x))
#log.debug("upper constr value {}".format(-x))
return -x
@ -178,7 +212,8 @@ def opt_integral_boundaries(integrand, a, b, t_max, ft_ref, opt_b_only, N, diff_
_f_opt_cache = dict()
args = (integrand, a, b, N, t_max, ft_ref, diff_method, _f_opt_cache, opt_b_only)
r = minimize(_f_opt, x0 = [-0.1], args = args,
x0 = np.log10(0.1*integrand(b))
r = minimize(_f_opt, x0 = x0, args = args,
method='SLSQP',
constraints=[{"type": "ineq", "fun": _lower_contrs, "args": args},
{"type": "ineq", "fun": _upper_contrs}])
@ -219,7 +254,7 @@ def get_N_a_b_for_accurate_fourier_integral(integrand, a, b, t_max, tol, ft_ref,
raise RuntimeError("maximum number of points for Fourier Transform reached")
i += 1
def get_dt_for_accurate_interpolation(t_max, tol, ft_ref):
def get_dt_for_accurate_interpolation(t_max, tol, ft_ref, diff_method=_absDiff):
N = 32
sub_sampl = 8
@ -231,43 +266,40 @@ def get_dt_for_accurate_interpolation(t_max, tol, ft_ref):
ft_intp = ComplexInterpolatedUnivariateSpline(x = tau_sub, y = ft_ref_n[::sub_sampl], k=3)
ft_intp_n = ft_intp(tau)
d = np.max(np.abs(ft_intp_n-ft_ref_n))
d = diff_method(ft_intp_n, ft_ref_n)
log.debug("acc interp N {} dt {:.3e} {:.3e} -> d {:.3e}".format(N, sub_sampl*tau[1], t_max/(N/sub_sampl), d))
if d < tol:
return t_max/(N/sub_sampl)
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):
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):
N, a, b = get_N_a_b_for_accurate_fourier_integral(integrand, a, b,
t_max = t_max,
tol = intgr_tol,
ft_ref = ft_ref,
opt_b_only=opt_b_only,
N_max = N_max)
N_max = N_max,
diff_method=diff_method)
dt_tol = get_dt_for_accurate_interpolation(t_max = t_max,
tol = intpl_tol,
ft_ref = ft_ref)
ft_ref = ft_ref,
diff_method=diff_method)
dx = (b-a)/N
dt = 2*np.pi/dx/N
if dt <= dt_tol:
log.info("dt criterion fulfilled")
log.debug("dt criterion fulfilled")
return a, b, N, dx, dt
else:
log.info("down scale dx and dt to match new power of 2 N")
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)
#assert scale <= 1
scale = 1
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
if opt_b_only:
b = a + b_minus_a
else:
@ -275,13 +307,4 @@ def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_
b += delta/2
a -= delta/2
rd, a, b = opt_integral_boundaries(integrand=integrand, a=a, b=b, t_max=t_max, ft_ref=ft_ref,
opt_b_only=opt_b_only, N=N)
log.debug("rd after final optimization:{:.3e}".format(rd))
return a, b, N, dx_new, dt_new

2
stocproc/method_kle.py
View file

@ -365,7 +365,7 @@ def auto_ng(corr, t_max, ngfac=2, meth=get_mid_point_weights_times, tol=1e-3, di
the interpolation error is maximal when beeing in between the reference points.
5) Now calculate the deviation :math:`\Delta(n)` for sequential n starting at n=0. Stop if
:math:`\Delta(n) < tol`. If the deviation does not drop below tol for all :math:`0 \leq n < ng-1` increase
ng as follows :math:`ng = 2*ng-1` and start over at 1). (This update schema for ng asured that ng is odd
ng as follows :math:`ng = 2*ng-1` and start over at 1). (This update scheme for ng asured that ng is odd
which is needed for the 'simpson' and 'fourpoint' integration weights)
.. note::

14
stocproc/stocproc.py
View file

@ -268,16 +268,16 @@ class StocProc_KLE(_absStocProc):
class StocProc_FFT(_absStocProc):
r"""Simulate Stochastic Process using FFT method
This method works only for auto correlations functions of the form
This method uses the relation of the auto correlation to the spectral density
.. math:: \alpha(\tau) = \int_{\omega_\mathrm{min}}^{\omega_\mathrm{max}} \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
.. math:: \alpha(\tau) = \int \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
where :math:`J(\omega)` is a real non negative function, usually called spectral density.
Then the integral can be approximated by a discrete integration schema:
Then the integral can be approximated by a discrete integration scheme:
.. math:: \alpha(\tau) \approx \sum_{k=0}^{N-1} w_k \frac{J(\omega_k)}{\pi} e^{-\mathrm{i} k \omega_k \tau}
where the :math:`\omega_k` depend on the integration schema.
where the :math:`\omega_k` depend on the integration scheme.
For a process defined by
@ -294,14 +294,14 @@ class StocProc_FFT(_absStocProc):
\approx & \alpha(t-s)
\end{align}
To calculate :math:`Z(t)` the schema of the Discrete Fourier Transform (DFT) can be applied as follows:
To calculate :math:`Z(t)` the Discrete Fourier Transform (DFT) can be applied as follows:
.. math:: Z(t_l) = e^{-\mathrm{i}\omega_\mathrm{min} t_l} \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k e^{-\mathrm{i} 2 \pi \frac{k l}{N} \frac{\Delta \omega \Delta t}{ 2 \pi} N}
Here :math:`\omega_k` has to take the form :math:`\omega_k = \omega_\mathrm{min} + k \Delta \omega` and
:math:`\Delta \omega = (\omega_\mathrm{max} - \omega_\mathrm{min}) / N-1` which limits
the itegration schemas to those with equidistant weights.
For the DFT schema to be applicable :math:`\Delta t` has to be chosen such that
the itegration schemes to those with equidistant weights.
For the DFT scheme to be applicable :math:`\Delta t` has to be chosen such that
.. math:: 1 = \frac{\Delta \omega \Delta t}{2 \pi} N

94
tests/test_method_fft.py
View file

@ -1,22 +1,22 @@
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 os
import pathlib
import scipy.integrate as sp_int
from scipy.special import gamma as gamma_func
import stocproc as sp
import sys
from stocproc import method_fft
try:
import matplotlib.pyplot as plt
except ImportError:
print("matplotlib not found -> any plotting will crash")
p = pathlib.PosixPath(os.path.abspath(__file__))
sys.path.insert(0, str(p.parent.parent))
def test_find_integral_boundary():
def f(x):
return np.exp(-(x)**2)
@ -41,14 +41,14 @@ def test_find_integral_boundary():
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)/f2(1) -tol) < 1e-14
assert abs(f2(b)/f2(-1)-tol) < 1e-14
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)/f2(1)-tol) < 1e-14, "diff {}".format(abs(f2(a)/f2(1)-tol))
assert abs(f2(b)/f2(-1)-tol) < 1e-14, "diff {}".format(abs(f2(b)/f2(-1)-tol))
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
@ -57,16 +57,30 @@ def test_find_integral_boundary():
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
f3_refval = f3(5)
assert abs(f3(a)/f3_refval-tol) < 1e-14
assert abs(f3(b)/f3_refval-tol) < 1e-14
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)/f3_refval-tol) < 1e-14, "diff {}".format(abs(f3(a)-tol))
assert abs(f3(b)/f3_refval-tol) < 1e-14, "diff {}".format(abs(f3(b)-tol))
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
@ -324,44 +338,44 @@ def test_calc_abN():
tol = 1e-3
tmax=40
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)
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)
print(dt, tau[1])
idx = np.where(tau <= tmax)
ft_ref_tau = bcf_ref(tau[idx])
rd = np.max(np.abs(ft_tau[idx] - ft_ref_tau) / np.abs(ft_ref_tau))
print("rd {:.3e} <= tol {:.3e}".format(rd, tol))
assert rd < tol
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)
ft_intp_n = ft_intp(tau_fine)
d = np.max(np.abs(ft_intp_n - ft_ref_n))
print("d {:.3e} <= tol {:.3e}".format(d, tol))
assert d < tol
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
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=True)
# 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_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()