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

doc work in progress

Browse source
This commit is contained in:
Richard Hartmann 2018-03-26 15:16:01 +02:00
parent fa54046cab
commit d483d110d4
9 changed files with 120 additions and 69 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
docs/source/api.rst
View file

@ -15,8 +15,8 @@ and the classes implementing the ``_abcStocProc`` interface.
StocProc_KLE StocProc_KLE
StocProc_TS StocProc_TS
Function on the mudule level Function on mudule level
---------------------------- ------------------------
.. autofunction:: stocproc.loggin_setup .. autofunction:: stocproc.loggin_setup
.. autofunction:: stocproc.version .. autofunction:: stocproc.version

14
docs/source/index.rst
View file

@ -2,14 +2,18 @@ StocProc
======== ========
The StocProc module is a Python3 module allowing to sample The StocProc module is a Python3 module allowing to sample
Gaussian stochastic processes which are wide-sense stationary, Gaussian stochastic processes :math:`z(t)` which are wide-sense stationary,
continuous in time and complex valued. In particular for a given auto correlation function continuous in time and complex valued. In particular for a given auto correlation function
:math:`\alpha(\tau)` the stochastic process :math:`z(t)` obeys. :math:`\alpha(\tau)` the stochastic process obeys:
.. math:: \langle z(t) \rangle = 0 \qquad \langle z(t)z(s) \rangle = 0 \qquad \langle z(t)z^\ast(s) \rangle = \alpha(t-s) .. math:: \langle z(t) \rangle = 0 \qquad \langle z(t)z(s) \rangle = 0 \qquad \langle z(t)z^\ast(s) \rangle = \alpha(t-s)
Here :math:`\langle \cdot \rangle` denotes the ensemble average. Here :math:`\langle \cdot \rangle` denotes the ensemble average.
The so far implemented methods (KLE, FFT, TanhSinh) to generate such processes provide an error control mechanism
which ensures that the auto correlation function of the numerically generated stochastic processes does correspond
to the preset auto correlation function :math:`\alpha(\tau)`.
Example Example
------- -------
@ -35,9 +39,9 @@ and sample a single realization.
print("setup process generator") print("setup process generator")
stp = sp.StocProc_FFT(spectral_density = lsd, stp = sp.StocProc_FFT(spectral_density = lsd,
t_max = t_max, t_max = t_max,
bcf_ref = exp_ac, alpha = exp_ac,
intgr_tol=1e-2, intgr_tol=1e-2, # integration error control parameter
intpl_tol=1e-2): intpl_tol=1e-2): # interpolation error control parameter
print("generate single process") print("generate single process")
stp.new_process() stp.new_process()

27
docs/source/methods.rst
View file

@ -35,6 +35,24 @@ For complex valued Gaussian distributed and independent random variables :math:`
obeys the required statistic. Note, the upper integral limit :math:`T` sets the time interval for which the obeys the required statistic. Note, the upper integral limit :math:`T` sets the time interval for which the
stochastic process :math:`z(t) \; t \in [0,T]` is defined. stochastic process :math:`z(t) \; t \in [0,T]` is defined.
In principal the sum is infinite. Nonetheless, a finite subset of summands can be found to yield a very good
approximation of the preset auto correlations functions.
Secondly when solving the Fredholm equation numerically, the integral is approximated in terms of a sum with
integration weights :math:`w_i`,
which in turn yields a matrix Eigenvalue problem with discrete "Eigenfunctions"
([NumericalRecipes]_ Chap. 19.1).
Comparing the preset auto correlation function with the approximate auto correlation function
using a finite set of :math:`N` discrete Eigenfunctions
.. math:: \sum_{n=1}^N \lambda_n u_n(t) u_n^\ast(s)
where :math:`u_n(t)` is the interpolated discrete Eigenfunction ([NumericalRecipes]_ eq. 19.1.3)
.. math:: u_n(t) = \sum_i \frac{w_i}{\lambda_n} \alpha(t-s_i) u_{n,i}
allows for an error estimation.
The KLE approach is implemented by the class :py:class:`stocproc.StocProc_KLE`. The KLE approach is implemented by the class :py:class:`stocproc.StocProc_KLE`.
It is numerically feasible if :math:`T` is not too large in comparision to a typical decay time of the It is numerically feasible if :math:`T` is not too large in comparision to a typical decay time of the
auto correlation function. auto correlation function.
@ -69,6 +87,12 @@ For complex valued Gaussian distributed and independent random variables :math:`
obeys the required statistics up to an accuracy of the integral discretization. obeys the required statistics up to an accuracy of the integral discretization.
To ensure efficient evaluation of the stochastic process the continuous time property is realized only approximately
by interpolating a pre calculated discrete time process.
However, the error caused by the cubic spline interpolation can be explicitly controlled
(usually by the `intpl_tol` parameter). Error values of one percent and below are easily achievable.
Fast Fourier Transform (FFT) Fast Fourier Transform (FFT)
```````````````````````````` ````````````````````````````
@ -80,3 +104,6 @@ TanhSinh Intgeration (TanhSinh)
For spectral densities :math:`J(\omega)` with a singularity at :math:`\omega=0` the TanhSinh integration For spectral densities :math:`J(\omega)` with a singularity at :math:`\omega=0` the TanhSinh integration
scheme is more suitable. Such an implementation and its details can be found at :py:class:`stocproc.StocProc_TanhSinh`. scheme is more suitable. Such an implementation and its details can be found at :py:class:`stocproc.StocProc_TanhSinh`.
.. [NumericalRecipes] Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes 3rd Edition: The Art of Scientific Computing, Auflage: 3. ed. Cambridge University Press, Cambridge, UK ; New York.

4
setup.py
View file

@ -1,13 +1,13 @@
from setuptools import setup from setuptools import setup
from Cython.Build import cythonize from Cython.Build import cythonize
import numpy as np import numpy as np
from stocproc import __version__ from stocproc import version_full
author = u"Richard Hartmann" author = u"Richard Hartmann"
authors = [author] authors = [author]
description = 'Generate continuous time stationary stochastic processes from a given auto correlation function.' description = 'Generate continuous time stationary stochastic processes from a given auto correlation function.'
name = 'stocproc' name = 'stocproc'
version = __version__ version = version_full()
if __name__ == "__main__": if __name__ == "__main__":
setup( setup(

3
stocproc/__init__.py
View file

@ -14,7 +14,10 @@ import sys
if sys.version_info.major < 3: if sys.version_info.major < 3:
raise SystemError("no support for Python 2") raise SystemError("no support for Python 2")
try:
from .stocproc import loggin_setup from .stocproc import loggin_setup
from .stocproc import StocProc_FFT from .stocproc import StocProc_FFT
from .stocproc import StocProc_KLE from .stocproc import StocProc_KLE
from .stocproc import StocProc_TanhSinh from .stocproc import StocProc_TanhSinh
except ImportError:
print("WARNING: Import Error occurred, parts of the package may be not available")

21
stocproc/method_ft.py
View file

@ -410,6 +410,27 @@ def get_dt_for_accurate_interpolation(t_max, tol, ft_ref, diff_method=_absDiff):
def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, ft_ref, opt_b_only, diff_method=_absDiff): def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, ft_ref, opt_b_only, diff_method=_absDiff):
r"""
Calculate the parameters for FFT method such that the error tolerance is met.
Free parameters are:
- :math:`\omega_{max}`: the upper bound of the Fourier integral. This parameter has to be large enough
In particular two criterion have to be met.
1)
:param integrand:
:param intgr_tol:
:param intpl_tol:
:param t_max:
:param ft_ref:
:param opt_b_only:
:param diff_method:
:return:
"""
log.info("get_dt_for_accurate_interpolation, please wait ...") log.info("get_dt_for_accurate_interpolation, please wait ...")
try: try:

8
stocproc/stocproc.py
View file

@ -217,7 +217,7 @@ class StocProc_KLE(_abcStocProc):
def __init__(self, alpha, t_max, tol=1e-2, ng_fac=4, meth='fourpoint', diff_method='full', dm_random_samples=10**4, def __init__(self, alpha, t_max, tol=1e-2, ng_fac=4, meth='fourpoint', diff_method='full', dm_random_samples=10**4,
seed=None, align_eig_vec=False, scale=1): seed=None, align_eig_vec=False, scale=1):
""" r"""
:param r_tau: the idesired auto correlation function of a single parameter tau :param r_tau: the idesired auto correlation function of a single parameter tau
:param t_max: specifies the time interval [0, t_max] for which the processes in generated :param t_max: specifies the time interval [0, t_max] for which the processes in generated
:param tol: maximal deviation of the auto correlation function of the sampled processes from :param tol: maximal deviation of the auto correlation function of the sampled processes from
@ -342,8 +342,8 @@ class StocProc_FFT(_abcStocProc):
This is ensured by automatically determining the number of sumands N and the integral This is ensured by automatically determining the number of sumands N and the integral
boundaries :math:`\omega_\mathrm{min}` and :math:`\omega_\mathrm{max}` such that boundaries :math:`\omega_\mathrm{min}` and :math:`\omega_\mathrm{max}` such that
discrete Fourier transform of the spectral density matches the desired auto correlation function discrete Fourier transform of the spectral density matches the preset auto correlation function
within the tolerance intgr_tol for all discrete :math:`t_l \in [0, t_{max}]`. within the tolerance `intgr_tol` for all discrete :math:`t_l \in [0, t_{max}]`.
As the time continuous process is generated via cubic spline interpolation, the deviation As the time continuous process is generated via cubic spline interpolation, the deviation
due to the interpolation is controlled by the parameter ``intpl_tol``. The maximum time step :math:`\Delta t` due to the interpolation is controlled by the parameter ``intpl_tol``. The maximum time step :math:`\Delta t`
@ -354,7 +354,7 @@ class StocProc_FFT(_abcStocProc):
criterion from the interpolation is met. criterion from the interpolation is met.
See :py:func:`stocproc.method_ft.calc_ab_N_dx_dt` for implementation details on how the See :py:func:`stocproc.method_ft.calc_ab_N_dx_dt` for implementation details on how the
tolerance criterion are met. Since the pre calculation may become time consuming the :py:class:`StocProc_FFT` tolerance criterion is met. Since the pre calculation may become time consuming the :py:class:`StocProc_FFT`
class can be pickled and unpickled. To identify a particular instance a unique key is formed by the tuple class can be pickled and unpickled. To identify a particular instance a unique key is formed by the tuple
``(alpha, t_max, intgr_tol, intpl_tol)``. ``(alpha, t_max, intgr_tol, intpl_tol)``.
It is advisable to use :py:func:`get_key` with keyword arguments to generate such a tuple. It is advisable to use :py:func:`get_key` with keyword arguments to generate such a tuple.

68
tests/test_method_fft.py
View file

@ -10,7 +10,7 @@ import numpy as np
import scipy.integrate as sp_int import scipy.integrate as sp_int
from scipy.special import gamma as gamma_func from scipy.special import gamma as gamma_func
import stocproc as sp import stocproc as sp
from stocproc import method_fft from stocproc import method_ft
try: try:
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -22,14 +22,14 @@ def test_find_integral_boundary():
return np.exp(-(x)**2) return np.exp(-(x)**2)
tol = 1e-10 tol = 1e-10
b = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=+1, max_val=1e6) b = sp.method_ft.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) a = sp.method_ft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=-1, max_val=1e6)
assert a != b assert a != b
assert abs(f(a)-tol) < 1e-14 assert abs(f(a)-tol) < 1e-14
assert abs(f(b)-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) b = sp.method_ft.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) a = sp.method_ft.find_integral_boundary(integrand=f, tol=tol, ref_val=0, x0=a-5, max_val=1e6)
assert a != b assert a != b
assert abs(f(a)-tol) < 1e-14 assert abs(f(a)-tol) < 1e-14
assert abs(f(b)-tol) < 1e-14 assert abs(f(b)-tol) < 1e-14
@ -38,14 +38,14 @@ def test_find_integral_boundary():
return np.exp(-(x)**2)*x**2 return np.exp(-(x)**2)*x**2
tol = 1e-10 tol = 1e-10
b = sp.method_fft.find_integral_boundary(integrand=f2, tol=tol, ref_val=1, x0=+1, max_val=1e6) b = sp.method_ft.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) a = sp.method_ft.find_integral_boundary(integrand=f2, tol=tol, ref_val=-1, x0=-1, max_val=1e6)
assert a != b assert a != b
assert abs(f2(a) -tol) < 1e-14 assert abs(f2(a) -tol) < 1e-14
assert abs(f2(b)-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) b = sp.method_ft.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) a = sp.method_ft.find_integral_boundary(integrand=f2, tol=tol, ref_val=-1, x0=a-5, max_val=1e6)
assert a != b assert a != b
assert abs(f2(a)-tol) < 1e-14, "diff {}".format(abs(f2(a)/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)) assert abs(f2(b)-tol) < 1e-14, "diff {}".format(abs(f2(b)/f2(-1)-tol))
@ -54,14 +54,14 @@ def test_find_integral_boundary():
return np.exp(-(x-5)**2)*x**2 return np.exp(-(x-5)**2)*x**2
tol = 1e-10 tol = 1e-10
b = sp.method_fft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=+1, max_val=1e6) b = sp.method_ft.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) a = sp.method_ft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=-1, max_val=1e6)
assert a != b assert a != b
assert abs(f3(a)-tol) < 1e-14 assert abs(f3(a)-tol) < 1e-14
assert abs(f3(b)-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) b = sp.method_ft.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) a = sp.method_ft.find_integral_boundary(integrand=f3, tol=tol, ref_val=5, x0=a-5, max_val=1e6)
assert a != b assert a != b
assert abs(f3(a)-tol) < 1e-14, "diff {}".format(abs(f3(a)-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)) assert abs(f3(b)-tol) < 1e-14, "diff {}".format(abs(f3(b)-tol))
@ -75,9 +75,9 @@ def test_find_integral_boundary():
return np.exp(-x ** 2) return np.exp(-x ** 2)
tol = 1e-3 tol = 1e-3
b = sp.method_fft.find_integral_boundary(integrand=f, tol=tol, ref_val=10, x0=+1, max_val=1e6) b = sp.method_ft.find_integral_boundary(integrand=f, tol=tol, ref_val=10, x0=+1, max_val=1e6)
assert abs(f(b) - tol) < 1e-14 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) a = sp.method_ft.find_integral_boundary(integrand=f, tol=tol, ref_val=-10, x0=-1., max_val=1e6)
assert abs(f(a) - tol) < 1e-14 assert abs(f(a) - tol) < 1e-14
@ -85,7 +85,7 @@ def fourier_integral_trapz(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
this function is here and not in method_fft because it has almost no this function is here and not in method_ft because it has almost no
advantage over the modpoint method. so only for testing purposes. advantage over the modpoint method. so only for testing purposes.
""" """
yl = integrand(np.linspace(a, b, N+1, endpoint=True)) yl = integrand(np.linspace(a, b, N+1, endpoint=True))
@ -151,7 +151,7 @@ def test_fourier_integral_finite_boundary():
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 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 = 2**18
N_test = 100 N_test = 100
tau, ft_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) tau, ft_n = sp.method_ft.fourier_integral_midpoint(intg, a, b, N)
ft_ref_n = ft_ref(tau) ft_ref_n = ft_ref(tau)
tau = tau[1:N_test] tau = tau[1:N_test]
ft_n = ft_n[1:N_test] ft_n = ft_n[1:N_test]
@ -174,7 +174,7 @@ def test_fourier_integral_finite_boundary():
## check against numeric fourier integral (non FFT) ## ## check against numeric fourier integral (non FFT) ##
###################################################### ######################################################
N = 512 N = 512
tau, ft_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) tau, ft_n = sp.method_ft.fourier_integral_midpoint(intg, a, b, N)
k, ft_simple = fourier_integral_simple_test(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 assert np.max(np.abs(ft_simple-ft_n)) < 1e-11
@ -189,7 +189,7 @@ def test_fourier_integral_finite_boundary():
## check midp against simpson ## ## check midp against simpson ##
################################# #################################
N = 1024 N = 1024
tau, ft_n = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) tau, ft_n = sp.method_ft.fourier_integral_midpoint(intg, a, b, N)
ft_ref_n = ft_ref(tau) ft_ref_n = ft_ref(tau)
rd = np.abs(ft_ref_n-ft_n) / np.abs(ft_ref_n) rd = np.abs(ft_ref_n-ft_n) / np.abs(ft_ref_n)
idx = np.where(np.logical_and(tau < 75, np.isfinite(rd))) idx = np.where(np.logical_and(tau < 75, np.isfinite(rd)))
@ -198,7 +198,7 @@ def test_fourier_integral_finite_boundary():
assert mrd_midp < 9e-3, "mrd_midp = {}".format(mrd_midp) assert mrd_midp < 9e-3, "mrd_midp = {}".format(mrd_midp)
N = 513 N = 513
tau, ft_n = sp.method_fft.fourier_integral_simps(intg, a, b, N) tau, ft_n = sp.method_ft.fourier_integral_simps(intg, a, b, N)
ft_ref_n = ft_ref(tau) ft_ref_n = ft_ref(tau)
rd = np.abs(ft_ref_n-ft_n) / np.abs(ft_ref_n) rd = np.abs(ft_ref_n-ft_n) / np.abs(ft_ref_n)
idx = np.where(np.logical_and(tau < 75, np.isfinite(rd))) idx = np.where(np.logical_and(tau < 75, np.isfinite(rd)))
@ -231,11 +231,11 @@ def test_fourier_integral_infinite_boundary(plot=False):
intg = lambda x: osd(x, s, wc) intg = lambda x: osd(x, s, wc)
bcf_ref = lambda t: gamma_func(s + 1) * wc**(s+1) * (1 + 1j*wc * t)**(-(s+1)) 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) a,b = sp.method_ft.find_integral_boundary_auto(integrand=intg, tol=1e-12, ref_val=1)
errs = [9e-5, 2e-5, 1.3e-6] errs = [9e-5, 2e-5, 1.3e-6]
for i, N in enumerate([2**16, 2**18, 2**20]): 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) tau, bcf_n = sp.method_ft.fourier_integral_midpoint(intg, a, b, N=N)
bcf_ref_n = bcf_ref(tau) bcf_ref_n = bcf_ref(tau)
tau_max = 5 tau_max = 5
@ -249,7 +249,7 @@ def test_fourier_integral_infinite_boundary(plot=False):
p, = plt.plot(tau, rd_mp, label="trapz N {}".format(N)) 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) tau, bcf_n = sp.method_ft.fourier_integral_simps(intg, a, b=b, N=N-1)
bcf_ref_n = bcf_ref(tau) bcf_ref_n = bcf_ref(tau)
idx = np.where(tau <= tau_max) idx = np.where(tau <= tau_max)
@ -283,8 +283,8 @@ def test_get_N_a_b_for_accurate_fourier_integral():
_WC_ = 2 _WC_ = 2
intg = lambda w: 1 / (1 + (w - _WC_) ** 2) / np.pi intg = lambda w: 1 / (1 + (w - _WC_) ** 2) / np.pi
bcf_ref = lambda t: np.exp(- np.abs(t) - 1j * _WC_ * t) 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_) a, b = sp.method_ft.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, N, a, b = sp.method_ft.get_N_a_b_for_accurate_fourier_integral(intg,
t_max=50, t_max=50,
tol=1e-2, tol=1e-2,
ft_ref=bcf_ref, ft_ref=bcf_ref,
@ -297,9 +297,9 @@ def test_get_N_a_b_for_accurate_fourier_integral_b_only():
intg = lambda x: osd(x, s, wc) intg = lambda x: osd(x, s, wc)
bcf_ref = lambda t: gamma_func(s + 1) * wc**(s+1) * (1 + 1j*wc * t)**(-(s+1)) 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, b = sp.method_ft.find_integral_boundary_auto(integrand=intg, tol=1e-2, ref_val=1)
a = 0 a = 0
N, a, b = sp.method_fft.get_N_a_b_for_accurate_fourier_integral(intg, N, a, b = sp.method_ft.get_N_a_b_for_accurate_fourier_integral(intg,
t_max=15, t_max=15,
tol=1e-5, tol=1e-5,
ft_ref=bcf_ref, ft_ref=bcf_ref,
@ -311,12 +311,12 @@ def test_get_dt_for_accurate_interpolation():
wc = 4 wc = 4
intg = lambda x: osd(x, s, wc) intg = lambda x: osd(x, s, wc)
bcf_ref = lambda t: gamma_func(s + 1) * wc**(s+1) * (1 + 1j*wc * t)**(-(s+1)) 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) dt = sp.method_ft.get_dt_for_accurate_interpolation(t_max=40, tol=1e-4, ft_ref=bcf_ref)
print(dt) print(dt)
def test_sclicing(): def test_sclicing():
yl = np.ones(10, dtype=int) yl = np.ones(10, dtype=int)
yl = sp.method_fft.get_fourier_integral_simps_weighted_values(yl) yl = sp.method_ft.get_fourier_integral_simps_weighted_values(yl)
assert yl[0] == 2/6 assert yl[0] == 2/6
assert yl[1] == 8/6 assert yl[1] == 8/6
assert yl[2] == 4/6 assert yl[2] == 4/6
@ -330,20 +330,16 @@ def test_sclicing():
def test_calc_abN(): def test_calc_abN():
def testing(intg, bcf_ref, tol, tmax): def testing(intg, bcf_ref, tol, tmax):
diff_method = method_fft._absDiff diff_method = method_ft._absDiff
a, b, N, dx, dt = sp.method_ft.calc_ab_N_dx_dt(integrand=intg,
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, intgr_tol=tol,
intpl_tol=tol, intpl_tol=tol,
t_max=tmax, t_max=tmax,
a=0,
b=b,
ft_ref=bcf_ref, ft_ref=bcf_ref,
opt_b_only=True, opt_b_only=True,
diff_method=diff_method) diff_method=diff_method)
tau, ft_tau = sp.method_fft.fourier_integral_midpoint(intg, a, b, N) tau, ft_tau = sp.method_ft.fourier_integral_midpoint(intg, a, b, N)
idx = np.where(tau <= tmax) idx = np.where(tau <= tmax)
ft_ref_tau = bcf_ref(tau[idx]) ft_ref_tau = bcf_ref(tau[idx])
rd = diff_method(ft_tau[idx], ft_ref_tau) rd = diff_method(ft_tau[idx], ft_ref_tau)

14
tests/test_stocproc.py
View file

@ -120,7 +120,7 @@ def test_stochastic_process_KLE_correlation_function(plot=False):
t_max = 15 t_max = 15
num_samples = 2000 num_samples = 2000
tol = 3e-2 tol = 3e-2
stp = sp.StocProc_KLE(tol=1e-2, r_tau=corr, t_max=t_max, ng_fac=4, seed=0) stp = sp.StocProc_KLE(tol=1e-2, alpha=corr, t_max=t_max, ng_fac=4, seed=0)
stocproc_metatest(stp, num_samples, tol, corr, plot) stocproc_metatest(stp, num_samples, tol, corr, plot)
@ -136,7 +136,7 @@ def test_stochastic_process_FFT_correlation_function(plot=False):
t_max = 15 t_max = 15
num_samples = 2000 num_samples = 2000
tol = 3e-2 tol = 3e-2
stp = sp.StocProc_FFT(spectral_density=spectral_density, t_max=t_max, bcf_ref=corr, intgr_tol=1e-2, intpl_tol=1e-2, stp = sp.StocProc_FFT(spectral_density=spectral_density, t_max=t_max, alpha=corr, intgr_tol=1e-2, intpl_tol=1e-2,
seed=0) seed=0)
stocproc_metatest(stp, num_samples, tol, corr, plot) stocproc_metatest(stp, num_samples, tol, corr, plot)
@ -147,7 +147,7 @@ def test_stocproc_dump_load():
## STOCPROC KLE ## ## STOCPROC KLE ##
#################### ####################
t0 = time.time() t0 = time.time()
stp = sp.StocProc_KLE(tol=1e-2, r_tau=corr, t_max=t_max, ng_fac=4, seed=0) stp = sp.StocProc_KLE(tol=1e-2, alpha=corr, t_max=t_max, ng_fac=4, seed=0)
t1 = time.time() t1 = time.time()
dt1 = t1 - t0 dt1 = t1 - t0
stp.new_process() stp.new_process()
@ -224,16 +224,16 @@ def test_many(plot=False):
stp = sp.StocProc_FFT(sd, t_max, ac, negative_frequencies=True, seed=0, intgr_tol=5e-3, intpl_tol=5e-3) stp = sp.StocProc_FFT(sd, t_max, ac, negative_frequencies=True, seed=0, intgr_tol=5e-3, intpl_tol=5e-3)
stocproc_metatest(stp, num_samples, tol, ac, plot) stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='simp') stp = sp.StocProc_KLE(tol=5e-3, alpha=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='simp')
stocproc_metatest(stp, num_samples, tol, ac, plot) stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='simp') stp = sp.StocProc_KLE(tol=5e-3, alpha=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='simp')
stocproc_metatest(stp, num_samples, tol, ac, plot) stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='fp') stp = sp.StocProc_KLE(tol=5e-3, alpha=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='fp')
stocproc_metatest(stp, num_samples, tol, ac, plot) stocproc_metatest(stp, num_samples, tol, ac, plot)
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='fp') stp = sp.StocProc_KLE(tol=5e-3, alpha=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='random', meth='fp')
stocproc_metatest(stp, num_samples, tol, ac, plot) stocproc_metatest(stp, num_samples, tol, ac, plot)
def test_pickle_scale(): def test_pickle_scale():