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

work in progress on cleanup and doc, getting ready for version 1.0.0

Browse source
This commit is contained in:
Richard Hartmann 2018-03-23 16:40:12 +01:00
parent 765c30636b
commit fa54046cab
15 changed files with 426 additions and 314 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
docs/source/StocProc_ABC.rst Normal file
View file

@ -0,0 +1,8 @@
Abstract Stochastic Process Class
---------------------------------
.. py:currentmodule:: stocproc.stocproc
.. autoclass:: stocproc.stocproc._abcStocProc
:members:
.. automethod:: __call__

5
docs/source/StocProc_FFT.rst
View file

@ -1,12 +1,9 @@
StocProc_FFT
------------
.. py:currentmodule:: stocproc.stocproc
.. autoclass:: stocproc.stocproc.StocProc_FFT
:members:
:inherited-members:
:undoc-members:
.. autoclass:: stocproc.stocproc.StocProc_TanhSinh
:members:
:inherited-members:
:undoc-members:

19
docs/source/StocProc_KLE.rst
View file

@ -10,12 +10,13 @@ StocProc_KLE
.. automethod:: __call__
.. py:currentmodule:: stocproc.method_kle
.. autofunction:: stocproc.method_kle.solve_hom_fredholm
.. autofunction:: stocproc.method_kle.get_mid_point_weights_times
.. autofunction:: stocproc.method_kle.get_trapezoidal_weights_times
.. autofunction:: stocproc.method_kle.get_simpson_weights_times
.. autofunction:: stocproc.method_kle.get_four_point_weights_times
.. autofunction:: stocproc.method_kle.get_gauss_legendre_weights_times
.. autofunction:: stocproc.method_kle.get_tanh_sinh_weights_times
.. autofunction:: stocproc.method_kle.auto_ng
..
.. py:currentmodule:: stocproc.method_kle
.. autofunction:: stocproc.method_kle.solve_hom_fredholm
.. autofunction:: stocproc.method_kle.get_mid_point_weights_times
.. autofunction:: stocproc.method_kle.get_trapezoidal_weights_times
.. autofunction:: stocproc.method_kle.get_simpson_weights_times
.. autofunction:: stocproc.method_kle.get_four_point_weights_times
.. autofunction:: stocproc.method_kle.get_gauss_legendre_weights_times
.. autofunction:: stocproc.method_kle.get_tanh_sinh_weights_times
.. autofunction:: stocproc.method_kle.auto_ng.

7
docs/source/StocProc_TS.rst Normal file
View file

@ -0,0 +1,7 @@
StocProc_TanhSinh
-----------------
.. autoclass:: stocproc.stocproc.StocProc_TanhSinh
:members:
:inherited-members:
:undoc-members:

34
docs/source/api.rst Normal file
View file

@ -0,0 +1,34 @@
API documentation
=================
The ``_abcStocProc`` interface definition
.. toctree::
StocProc_ABC
and the classes implementing the ``_abcStocProc`` interface.
.. toctree::
:maxdepth: 3
StocProc_FFT
StocProc_KLE
StocProc_TS
Function on the mudule level
----------------------------
.. autofunction:: stocproc.loggin_setup
.. autofunction:: stocproc.version
.. autofunction:: stocproc.version_full
Method related submodules
-------------------------
.. toctree::
:maxdepth: 3
method_kle
method_ft

7
docs/source/conf.py
View file

@ -68,9 +68,12 @@ author = 'Richard Hartmann'
# built documents.
#
# The short X.Y version.
version = '1.0'
import stocproc
version = stocproc.version()
# The full version, including alpha/beta/rc tags.
release = '1.0.0'
release = stocproc.version_full()
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

4
docs/source/example.rst
View file

@ -1,5 +1,5 @@
Example
=======
Full Example
============
This example sets up a generator for stochastic processes with exponential
auto correlation. A single process is shown as well as the reconstruction

78
docs/source/index.rst
View file

@ -1,2 +1,78 @@
.. automodule:: stocproc
StocProc
========
The StocProc module is a Python3 module allowing to sample
Gaussian stochastic processes which are wide-sense stationary,
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:: \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.
Example
-------
The example will setup a process generator for an exponential auto correlation function
and sample a single realization.
::
def lsd(w):
# Lorenzian spectral density
return 1/(1 + (w - _WC_)**2)
def exp_ac(t):
# exponential auto correlation function
# note there is a factor of one over pi in the definition of the auto correlation function
# exp_ac(t) = 1/pi int_{-infty}^infty d w lsd(w) exp(-i w t)
return np.exp(- np.abs(t) - 1j*_WC_*t)
_WC_ = 5
t_max = 10
print("setup process generator")
stp = sp.StocProc_FFT(spectral_density = lsd,
t_max = t_max,
bcf_ref = exp_ac,
intgr_tol=1e-2,
intpl_tol=1e-2):
print("generate single process")
stp.new_process()
zt = stp() # get discrete process
.. image:: ../../examples/proc.*
:width: 500px
:align: center
:height: 300px
The full example code can be found :doc:`here </example>`.
How to Install
--------------
Install the latest version via pip ::
pip install stocproc
or fetch the bleeding edge version from the `git repository <https://github.com/cimatosa/stocproc>`_ ::
git clone https://github.com/cimatosa/stocproc.git
and install the package invoking ``setup.py``.
Note: The stocproc module depends on `fcSpline <https://github.com/cimatosa/fcSpline>`_ , a fast cubic spline interpolator for equally distributed nodes.
Documentation
-------------
.. toctree::
:maxdepth: 3
methods
api
example

10
docs/source/method_ft.rst Normal file
View file

@ -0,0 +1,10 @@
method_ft Module
-----------------
helper functions for the FT method
.. py:currentmodule:: stocproc
.. automodule:: stocproc.method_ft
:members:
:undoc-members:

9
docs/source/method_kle.rst Normal file
View file

@ -0,0 +1,9 @@
method_kle Module
-----------------
helper functions for the KLE method
.. automodule:: stocproc.method_kle
:members:
:undoc-members:

82
docs/source/methods.rst Normal file
View file

@ -0,0 +1,82 @@
Methods to Sample Stochastic Processes
======================================
This module allows to sample Gaussian stochastic processes which are wide-sense stationary,
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:: \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.
So far two different approaches to sample such processes have been implemented.
Karhunen-Loève Expansion (KLE)
------------------------------
The Karhunen-Loève expansion makes use of the spectral representation of the auto correlation function viewed as a
semi-positive integral kernel which allows for the decomposition.
.. math:: \alpha(t-s) = \sum_n \lambda_n u_n(t) u_n^\ast(s)
The Eigenvalues :math:`\lambda_n` and the Eigenfunction :math:`u_n(t)` of :math:`\alpha(t-s)` are defined in
terms of the homogeneous Fredholm equation.
.. math:: \int_0^T \mathrm{d} s \; \alpha(t-s) u_n(s) = \lambda_n u_n(t)
For complex valued Gaussian distributed and independent random variables :math:`Y_n` with
:math:`\langle Y_n \rangle = 0`, :math:`\langle Y_n Y_m \rangle = 0` and
:math:`\langle Y_n Y^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
.. math:: z(t) = \sum_n Y_n \sqrt{\lambda_n} u_n(t)
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.
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
auto correlation function.
Implementation details can be found in the class documentation of :py:class:`stocproc.StocProc_KLE`.
::
import stocproc as sp
import numpy as np
stp = sp.StocProc_KLE(alpha = lambda t: np.exp(- np.abs(t) - 1j*5*t),
t_max = 2)
stp.new_process()
zt = stp()
Fourier Transform Methods
-------------------------
This approach utilizes the relation of the auto correlation function and its Fourier transform.
.. math:: \alpha(\tau) = \frac{1}{\pi} \int_{-\infty}^{\infty} \mathrm{d}\omega \; J(\omega) e^{-\mathrm{i}\omega\tau} \qquad J(\omega) \geq 0
Discretizing the integral yields an approximate expression for the auto correlation function.
.. math:: \alpha(\tau) \approx \sum_n w_n \frac{J(\omega_n)}{\pi} e^{-\mathrm{i}\omega_n\tau}
For complex valued Gaussian distributed and independent random variables :math:`Y_n` with
:math:`\langle Y_n \rangle = 0`, :math:`\langle Y_n Y_m \rangle = 0` and
:math:`\langle Y_n Y^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
.. math:: z(t) = \sum_n Y_n \sqrt{\frac{w_n J(\omega_n)}{\pi}} e^{-\mathrm{i}\omega_n t}
obeys the required statistics up to an accuracy of the integral discretization.
Fast Fourier Transform (FFT)
````````````````````````````
Equally distributed nodes :math:`\omega_n` allow for an evaluation of the stochastic process
using the Fast Fourier Transform algorithm (see :py:class:`stocproc.StocProc_FFT` for an implementation).
TanhSinh Intgeration (TanhSinh)
```````````````````````````````
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`.

1
docs/source/stocproc.rst
View file

@ -1 +0,0 @@
.. automodule:: stocproc.stocproc

150
stocproc/__init__.py
View file

@ -1,150 +1,20 @@
r"""
Stochastic Process Module
=========================
__MAJOR__ = 1
__MINOR__ = 0
__PATCH__ = 0
def version():
"""semantic version string with format 'MAJOR.MINOR' (https://semver.org/)"""
return "{}.{}".format(__MAJOR__, __MINOR__)
This module allows to sample Gaussian stochastic processes which are wide-sense stationary,
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:: \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.
So far two different approaches to sample such processes have been implemented.
:doc:`Karhunen-Loève expansion (KLE) </StocProc_KLE>`
-----------------------------------------------------
The Karhunen-Loève expansion makes use of the spectral representation of the auto correlation function viewed as a
semi-positive integral kernel which allows for the decomposition.
.. math:: \alpha(t-s) = \sum_n \lambda_n u_n(t) u_n^\ast(s)
The Eigenvalues :math:`\lambda_n` and the Eigenfunction :math:`u_n(t)` of :math:`\alpha(t-s)` are defined in
terms of the homogeneous Fredholm equation.
.. math:: \int_0^T \mathrm{d} s \; \alpha(t-s) u_n(s) = \lambda_n u_n(t)
For complex valued Gaussian distributed and independent random variables :math:`Z_n` with
:math:`\langle Z_n \rangle = 0`, :math:`\langle Z_n Z_m \rangle = 0` and
:math:`\langle Z_n Z^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
.. math:: z(t) = \sum_n Z_n \sqrt{\lambda_n} u_n(t)
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.
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
auto correlation function.
Implementation details can be found in the class documentation of :py:class:`stocproc.StocProc_KLE`.
::
import stocproc as sp
import numpy as np
stp = sp.StocProc_KLE(alpha = lambda t: np.exp(- np.abs(t) - 1j*5*t),
t_max = 2)
stp.new_process()
zt = stp()
:doc:`Fourier transform (FT) </StocProc_FFT>`
---------------------------------------------
This approach utilizes the relation of the auto correlation function and its Fourier transform.
.. math:: \alpha(\tau) = \frac{1}{\pi} \int_{-\infty}^{\infty} \mathrm{d}\omega \; J(\omega) e^{-\mathrm{i}\omega\tau} \qquad J(\omega) \geq 0
Discretizing the integral yields an approximate expression for the auto correlation function.
.. math:: \alpha(\tau) \approx \sum_n w_n \frac{J(\omega_n)}{\pi} e^{-\mathrm{i}\omega_n\tau}
For complex valued Gaussian distributed and independent random variables :math:`Z_n` with
:math:`\langle Z_n \rangle = 0`, :math:`\langle Z_n Z_m \rangle = 0` and
:math:`\langle Z_n Z^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
.. math:: z(t) = \sum_n Z_n \sqrt{\frac{w_n J(\omega_n)}{\pi}} e^{-\mathrm{i}\omega_n t}
obeys the required statistics up to an accuracy of the integral discretization.
Equally distributed nodes :math:`\omega_n` allow for an evaluation of the stochastic process
using the Fast Fourier Transform algorithm (see :py:class:`stocproc.StocProc_FFT` for an implementation).
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`.
To implement additional methods
This module contains two different implementation for generating stochastic processes for a
given auto correlation function (:doc:`Karhunen-Loève expansion </StocProc_KLE>`
and :doc:`Fast-Fourier method </StocProc_FFT>`).
Both methods are based on a time discrete process, however cubic
spline interpolation is assured to be valid within a given tolerance.
Documentation Overview
----------------------
.. toctree::
:maxdepth: 3
stocproc
example
a Simple Example
----------------
The example will setup a process generator for an exponential auto correlation function
and sample a single realization. ::
def lsd(w):
# Lorenzian spectral density
return 1/(1 + (w - _WC_)**2)
def exp_ac(t):
# exponential auto correlation function
# note there is a factor of one over pi in the
# definition of the auto correlation function:
# exp_ac(t) = 1/pi int_{-infty}^infty d w lsd(w) exp(-i w t)
return np.exp(- np.abs(t) - 1j*_WC_*t)
_WC_ = 5
t_max = 10
print("setup process generator")
stp = sp.StocProc_FFT(spectral_density = lsd,
t_max = t_max,
bcf_ref = exp_ac,
intgr_tol=1e-2,
intpl_tol=1e-2):
print("generate single process")
stp.new_process()
zt = stp() # get discrete process
The full example can be found :doc:`here </example>`.
.. image:: ../../examples/proc.*
Averaging over 5000 samples yields the auto correlation function which is in good agreement
with the exact auto correlation.
.. image:: ../../examples/ac.*
"""
__version__ = "0.2.1"
def version_full():
"""semantic version string with format 'MAJOR.MINOR.PATCH' (https://semver.org/)"""
return "{}.{}.{}".format(__MAJOR__, __MINOR__, __PATCH__)
import sys
if sys.version_info.major < 3:
raise SystemError("no support for Python 2")
from .stocproc import loggin_setup
from .stocproc import StocProc_FFT
from .stocproc import StocProc_KLE
from .stocproc import StocProc_TanhSinh

0
stocproc/method_fft.py → stocproc/method_ft.py
View file

326
stocproc/stocproc.py
View file

@ -1,93 +1,86 @@
"""
Stochastic Process Generators
=============================
Karhunen-Loève expansion
------------------------
This method samples stochastic processes using Karhunen-Loève expansion and
is implemented in the class :doc:`StocProc_KLE </StocProc_KLE>`.
Setting up the class involves solving an eigenvalue problem which grows with
the time interval the process is simulated on. Further generating a new process
involves a multiplication with that matrix, therefore it scales quadratically with the
time interval. Nonetheless it turns out that this method requires less random numbers
than the Fast-Fourier method.
Fast-Fourier method
-------------------
In the class :doc:`StocProc_FFT </StocProc_FFT>` a method based on Fast-Fourier transform is
used to sample stochastic processes.
Setting up this class is quite efficient as it only calculates values of the
associated spectral density. The number scales linear with the time interval of interest. However to achieve
sufficient accuracy many of these values are required. As the generation of a new process is based on
a Fast-Fouried-Transform over these values, this part is comparably lengthy.
Abstract Base Class for a Stochastic Process Generator
------------------------------------------------------
:doc:`_absStocProc </StocProc_ABC>`.
"""
import abc
import numpy as np
import time
from . import method_kle
from . import method_fft
from . import method_ft
import fcSpline
import logging
log = logging.getLogger(__name__)
ch = logging.StreamHandler()
ch.setFormatter(logging.Formatter('%(name)s - %(levelname)s - %(message)s'))
sh = logging.StreamHandler()
sh.setFormatter(logging.Formatter('%(name)s - %(levelname)s - %(message)s'))
def show_log(level = logging.INFO):
ch.setLevel(level)
def loggin_setup(sh_level = logging.INFO, sp_log_level = logging.INFO,
kle_log_level = logging.INFO, ft_log_level = logging.INFO):
"""
controls the logging levels
log.addHandler(ch)
log.setLevel(level)
On loading the ``stocproc`` module this function is called with its default argument
which ammount to informative logging.
method_kle.log.addHandler(ch)
method_kle.log.setLevel(level)
:param sh_level: the logging level for the (output) StreamHandler
:param sp_log_level: the logging level of the stocproc module log
:param kle_log_level: the logging level of the KLE helper functions log (see :doc:`method_kle`)
:param ft_log_level: the logging level of the FT helper functions log (see :doc:`method_ft`)
"""
sh.setLevel(sh_level)
method_fft.log.addHandler(ch)
method_fft.log.setLevel(level)
log.addHandler(sh)
log.setLevel(sp_log_level)
method_kle.log.addHandler(sh)
method_kle.log.setLevel(kle_log_level)
method_ft.log.addHandler(sh)
method_ft.log.setLevel(ft_log_level)
loggin_setup()
class _abcStocProc(abc.ABC):
r"""
Abstract base class to stochastic process interface
general work flow:
- Specify the time axis of interest [0, t_max] and it resolution (number of grid points), :math:`t_i = i \frac{t_max}{N_t-1}`.
- To evaluate the stochastic process at these points, a mapping from :math:`N_z` normal distributed
random complex numbers with :math:`\langle y_i y_j^\ast \rangle = 2 \delta_{ij}`
to the stochastic process :math:`z_{t_i}` is needed and depends on the implemented method (:py:func:`calc_z`).
- A new process should be generated by calling :py:func:`new_process`.
- When the __call__ method is invoked the results will be interpolated between the :math:`z_{t_i}`.
"""
def __init__(self, t_max=None, num_grid_points=None, seed=None, t_axis=None, scale=1):
r"""
:param t_max: specify time axis as [0, t_max], if None, the times must be explicitly
given by t_axis
:param num_grid_points: number of equidistant times on that axis
:param seed: if not ``None`` set seed to ``seed``
:param t_axis: an explicit definition of times t_k (may be non equidistant)
"""
if t_max is not None:
self.t = np.linspace(0, t_max, num_grid_points)
else:
self.t = t_axis
Interface definition for stochastic process implementations
self.t_max = self.t[-1]
self.num_grid_points = len(self.t)
A new implementation for a stochastic process generator should subclass :py:class:`_abcStocProc` and
overwrite :py:func:`calc_z` and :py:func:`get_num_y`.
Depending on the equally spaced times :math:`t_n = n \frac{t_{max}}{N-1}` with :math:`n = 0 \dots N-1`
(:math:`N` = number of grid points),
the function :py:func:`calc_z` should map :math:`M` independent complex valued and Gaussian
distributed random variables :math:`Y_m` (with :math:`\langle Y_m \rangle = 0 = \langle Y_m Y_{m'} \rangle` and
:math:`\langle Y_m Y^\ast_{m'} = \delta_{m m'}\rangle`) to the discrete time stochastic process :math:`z_n = z(t_n)`.
:py:func:`get_num_y` needs to return the number :math:`M` of random variables :math:`Y_m` required as
input for :py:func:`calc_z`.
Having implemented :py:func:`calc_z` and :py:func:`get_num_y` the :py:class:`_abcStocProc` provides
convenient functions such as:
- :py:func:`__call__`: evaluate the stochastic process for any time within the interval
:math:`[0, t_{max}]` using cubic spline interpolation
- :py:func:`get_time`: returns the times :math:`t_n`
- :py:func:`get_z`: returns the discrete stochastic process :math:`z_n`
- :py:func:`new_process`: draw new samples :math:`Y_m` and update :math:`z_n` as well as the
cubic spline interpolator
- :py:func:`set_scale`: set a scalar pre factor for the auto correlation function
which scales the stochastic process such that :math:`\langle z(t) z^\ast(s)\rangle = \text{scale} \; \alpha(t-s)`.
:param t_max: specifies the upper bound of the time interval
:param num_grid_points: number of grid points :math:`N`
:param seed: if not ``None`` seed the random number generator with ``seed``
:param t_axis: an explicit definition of times t_k (may be non equidistant)
:param scale: passes ``scale`` to :py:func:`set_scale`
Note: :py:func:`new_process` is **not** called on init. If you want to evaluate a particular
realization of the stocastic process, a new sample needs to be drawn by calling :py:func:`new_process`.
Otherwise a ``RuntimeError`` is raised.
"""
def __init__(self, t_max=None, num_grid_points=None, seed=None, scale=1):
self.t_max = t_max
self.num_grid_points = num_grid_points
self.t = np.linspace(0, t_max, num_grid_points)
self._z = None
self._interpolator = None
@ -101,13 +94,12 @@ class _abcStocProc(abc.ABC):
log.debug("init StocProc with t_max {} and {} grid points".format(t_max, num_grid_points))
def __call__(self, t=None):
r"""evaluates the stochastic process via spline interpolation of the discrete process :math:`z_k`
r"""Evaluates the stochastic process.
:param t: time to evaluate the stochastic process at, float of array of floats, if t is None
return the discrete process :math:`z_k` which corresponds to the times :math:`t_k` given by the
integration weights method
:return: a single complex value or a complex array of the shape of t that holds the values of
stochastic process at times t
If ``t`` is not ``None`` cubic spline interpolation is used to evaluate :math:`z(t)` based on the
discrete realization of the process :math:`z_n = z(t_n)`. ``t`` may be a single value or array like.
If ``t`` is ``None`` the discrete process :math:`z_n` is returned.
"""
if self._z is None:
raise RuntimeError("StocProc has NO random data, call 'new_process' to generate a new random process")
@ -119,10 +111,13 @@ class _abcStocProc(abc.ABC):
@abc.abstractmethod
def calc_z(self, y):
r"""
maps the normal distributed complex valued random variables y to the stochastic process
:return: the stochastic process, array of complex numbers
r"""*abstract method*
An implementation needs to map :math:`M` independent complex valued and Gaussian
distributed random variables :math:`Y_m` (with :math:`\langle Y_m \rangle = 0 = \langle Y_m Y_{m'} \rangle` and
:math:`\langle Y_m Y^\ast_{m'} = \delta_{m m'}\rangle`) to the discrete time stochastic process :math:`z_n = z(t_n)`.
:return: the discrete time stochastic process :math:`z_n`, array of complex numbers
"""
pass
@ -132,41 +127,46 @@ class _abcStocProc(abc.ABC):
@abc.abstractmethod
def get_num_y(self):
r"""
:return: number of complex random variables needed to calculate the stochastic process
r"""*abstract method*
An implementation needs to return the number :math:`M` of random variables :math:`Y_m` required as
input for :py:func:`calc_z`.
"""
pass
def get_time(self):
r"""Returns the time :math:`t_k` corresponding to the values :math:`z_k`
These times are determined by the integration weights method.
"""
r"""Returns the times :math:`t_n` for which the discrete time stochastic process may be evaluated."""
return self.t
def get_z(self):
r"""Returns the discrete process :math:`z_k`."""
r"""Returns the discrete time stochastic process :math:`z_n = z(t_n)`."""
return self._z
def new_process(self, y=None, seed=None):
r"""generate a new process by evaluating :py:func:`calc_z` with new random variables :math:`Y_i`
r"""Generate a new realization of the stochastic process.
:param y: independent normal distributed complex valued random variables with :math:`\sigma_{ij}^2 = \langle y_i y_j^\ast \rangle = 2 \delta_{ij}`
:param seed: if not None set seed to seed before generating samples
When y is given use these random numbers as input for :py:func:`calc_z`
otherwise generate a new set of random numbers.
If ``seed`` is not ``None`` seed the random number generator with ``seed`` before drawing new random numbers.
If ``y`` is ``None`` draw new random numbers to generate the new realization. Otherwise use ``y`` as input to
generate the new realization.
"""
t0 = time.time()
self._interpolator = None
# clean up old data
del self._interpolator
del self._z
self._proc_cnt += 1
if seed != None:
log.info("use fixed seed ({})for new process".format(seed))
log.info("use fixed seed ({}) for new process".format(seed))
np.random.seed(seed)
if y is None:
#random complex normal samples
y = np.random.normal(scale=self._one_over_sqrt_2, size = 2*self.get_num_y()).view(np.complex)
del self._z
else:
if len(y) != self.get_num_y():
raise RuntimeError("the length of 'y' ({}) needs to be {}".format(len(y), self.get_num_y()))
self._z = self._calc_scaled_z(y)
log.debug("proc_cnt:{} new process generated [{:.2e}s]".format(self._proc_cnt, time.time() - t0))
t0 = time.time()
@ -174,6 +174,10 @@ class _abcStocProc(abc.ABC):
log.debug("created interpolator [{:.2e}s]".format(time.time() - t0))
def set_scale(self, scale):
r"""
Set a scalar pre factor for the auto correlation function
which scales the stochastic process such that :math:`\langle z(t) z^\ast(s)\rangle = \text{scale} \; \alpha(t-s)`.
"""
self.scale = scale
self.sqrt_scale = np.sqrt(scale)
@ -294,18 +298,19 @@ class StocProc_KLE(_abcStocProc):
class StocProc_FFT(_abcStocProc):
r"""Simulate Stochastic Process using FFT method
r"""Generate Stochastic Processes using the Fast Fourier Transform (FFT) method
This method uses the relation of the auto correlation to the non negative real valued
spectral density :math:`J(\omega)`. The integral can be approximated by a discrete integration scheme
This method uses the relation of the auto correlation function ``alpha`` to the non negative real valued
spectral density (``spectral_density``) :math:`J(\omega)`.
The integral can be approximated by a discrete integration scheme
.. math::
\alpha(\tau) = \int_{\omega_\mathrm{min}}^{\omega_\mathrm{max}} \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
\alpha(\tau) = \int_{-\infty}^{\infty} \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
\approx \sum_{k=0}^{N-1} w_k \frac{J(\omega_k)}{\pi} e^{-\mathrm{i} \omega_k \tau}
where the weights :math:`\omega_k` depend on the particular integration scheme. For a process defined as
.. math:: Z(t) = \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \exp^{-\mathrm{i}\omega_k t}
.. math:: z(t) = \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \exp^{-\mathrm{i}\omega_k t}
with independent complex random variables :math:`Y_k` such that :math:`\langle Y_k \rangle = 0`,
:math:`\langle Y_k Y_{k'}\rangle = 0` and :math:`\langle Y_k Y^\ast_{k'}\rangle = \delta_{k,k'}`
@ -313,90 +318,97 @@ class StocProc_FFT(_abcStocProc):
.. math::
\begin{align}
\langle Z(t) Z^\ast(s) \rangle = & \sum_{k,k'} \frac{1}{\pi} \sqrt{w_k w_{k'} J(\omega_k)J(\omega_{k'})} \langle Y_k Y_{k'}\rangle \exp(-\mathrm{i}(\omega_k t - \omega_k' s)) \\
= & \sum_{k} \frac{w_k}{\pi} J(\omega_k) e^{-\mathrm{i}\omega_k (t-s)}
\langle z(t) z^\ast(s) \rangle = & \sum_{k,k'} \frac{1}{\pi} \sqrt{w_k w_{k'} J(\omega_k)J(\omega_{k'})} \langle Y_k Y_{k'}\rangle \exp(-\mathrm{i}(\omega_k t - \omega_k' s)) \\
= & \sum_{k} \frac{w_k}{\pi} J(\omega_k) e^{-\mathrm{i}\omega_k (t-s)} \\
\approx & \alpha(t-s)
\end{align}
To calculate :math:`Z(t)` 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}
.. 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}
However, this requires that :math:`\omega_k` takes the form :math:`\omega_k = \omega_\mathrm{min} + k \Delta \omega`
with :math:`\Delta \omega = (\omega_\mathrm{max} - \omega_\mathrm{min}) / (N-1)` which limits
the integration schemes to those with equidistant nodes.
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 schemes to those with equidistant weights.
For the DFT scheme to be applicable :math:`\Delta t` has to be chosen such that
:math:`2\pi = N \Delta \omega \Delta t` holds.
Since :math:`J(\omega)` is real it follows that :math:`X(t_l) = X^\ast(t_{N-l})`.
Since :math:`J(\omega)` is real it follows that :math:`z(t_l) = z^\ast(t_{N-l})`.
For that reason the stochastic process has only :math:`(N+1)/2` (odd :math:`N`) and
:math:`(N/2 + 1)` (even :math:`N`) independent time grid points.
To generate a process with given auto correlation function on the interval [0, t_max]
requires that the auto correlation function approximation is valid for all t in [0, t_max].
To generate a process with given auto correlation function on the interval :math:`[0, t_{max}]`
requires that the auto correlation function approximation is valid for all :math:`t` in :math:`[0, t_{max}]`.
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
discrete Fourier transform of the spectral density matches the desired auto correlation function
within the tolerance intgr_tol for all discrete :math:`t_l \in [0, t_\mathrm{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
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`
is chosen such that the interpolated valued at each half step :math:`t_i + \Delta t /2` differs at
most intpl_tol from the exact value of the auto correlation function.
most ``intpl_tol`` from the exact value of the auto correlation function.
If not fulfilled already N and the integration boundaries are increased such that the :math:`\Delta t`
criterion from the interpolation is met.
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`
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)``.
It is advisable to use :py:func:`get_key` with keyword arguments to generate such a tuple.
:param spectral_density: the spectral density :math:`J(\omega)` as callable function object
:param t_max: :math:`[0,t_\mathrm{max}]` is the interval for which the process will be calculated
:param bcf_ref: a callable which evaluates the Fourier integral exactly
:param alpha: a callable which evaluates the Fourier integral exactly
:param intgr_tol: tolerance for the integral approximation
:param intpl_tol: tolerance for the interpolation
:param seed: if not None, use this seed to seed the random number generator
:param negative_frequencies: if False, keep :math:`\omega_\mathrm{min} = 0` otherwise
find a negative :math:`\omega_\mathrm{min}` appropriately just like :math:`\omega_\mathrm{max}`
.. todo::
implement bcf_ref = None and use numeric integration as default
"""
def __init__(self, spectral_density, t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2,
def __init__(self, spectral_density, t_max, alpha, 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
self.key = self.get_key(alpha=alpha, t_max=t_max, intgr_tol=intgr_tol, intpl_tol=intpl_tol)
if not negative_frequencies:
log.info("non neg freq only")
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,
ft_ref = lambda tau:bcf_ref(tau)*np.pi,
opt_b_only= True)
a, b, N, dx, dt = method_ft.calc_ab_N_dx_dt(integrand = spectral_density,
intgr_tol = intgr_tol,
intpl_tol = intpl_tol,
t_max = t_max,
ft_ref = lambda tau:alpha(tau)*np.pi,
opt_b_only= True)
else:
log.info("use neg freq")
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,
ft_ref = lambda tau:bcf_ref(tau)*np.pi,
opt_b_only= False)
a, b, N, dx, dt = method_ft.calc_ab_N_dx_dt(integrand = spectral_density,
intgr_tol = intgr_tol,
intpl_tol = intpl_tol,
t_max = t_max,
ft_ref = lambda tau:alpha(tau)*np.pi,
opt_b_only= False)
assert abs(2*np.pi - N*dx*dt) < 1e-12
d = abs(2*np.pi - N*dx*dt)
if d >= 1e-12:
log.fatal('method_ft.calc_ab_N_dx_dt returned inconsistent data!')
raise RuntimeError('d = {:.3e} < 1e-12 FAILED!'.format(d))
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))
log.info("Fourier Integral Boundaries: [{:.3e}, {:.3e}]".format(a,b))
log.info("Number of Nodes : {}".format(N))
log.info("yields dx : {:.3e}".format(dx))
log.info("yields dt : {:.3e}".format(dt))
log.info("yields t_max : {:.3e}".format( (N-1)*dt))
num_grid_points = int(np.ceil(t_max/dt))+1
assert num_grid_points <= N
if num_grid_points > N:
log.fatal("num_grid_points and number of points used for FFT inconsistent!")
raise RuntimeError("num_grid_points = {} <= N_DFT = {} FAILED!".format(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,
@ -407,8 +419,12 @@ class StocProc_FFT(_abcStocProc):
self.omega_min_correction = np.exp(-1j*(a+dx/2)*self.t) #self.t is from the parent class
@staticmethod
def get_key(t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2):
return bcf_ref, t_max, intgr_tol, intpl_tol
def get_key(t_max, alpha, intgr_tol=1e-2, intpl_tol=1e-2):
"""
Returns the tuple ``(alpha, t_max, intgr_tol, intpl_tol)`` which uniquely identifies a particular
:py:class:`StocProc_FFT` instance
"""
return alpha, t_max, intgr_tol, intpl_tol
def __getstate__(self):
return self.yl, self.num_grid_points, self.omega_min_correction, self.t_max, self._seed, self.scale, self.key
@ -421,20 +437,20 @@ class StocProc_FFT(_abcStocProc):
scale = scale)
def calc_z(self, y):
r"""calculate
r"""Calculate the discrete time stochastic process using FFT algorithm
.. math::
Z(t_l) = e^{-\mathrm{i}\omega_\mathrm{min} t_l} \mathrm{DFT}\left( \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \right)
z_n = z(t_n) = e^{-\mathrm{i}\omega_\mathrm{min} t_n} \mathrm{FFT}\left( \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \right)
and return values with :math:`t_l < t_\mathrm{max}`
and return values :math:`z_n` with :math:`t_n <= t_\mathrm{max}`.
"""
z_fft = np.fft.fft(self.yl * y)
z = z_fft[0:self.num_grid_points] * self.omega_min_correction
return z
def get_num_y(self):
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
r"""The number of independent random variables :math:`Y_m` is given by the number of discrete nodes
used by the Fast Fourier Transform algorithm.
"""
return len(self.yl)
@ -452,13 +468,13 @@ class StocProc_TanhSinh(_abcStocProc):
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)),
c = method_ft.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,
dt_tol = method_ft.get_dt_for_accurate_interpolation(t_max=c,
tol=intpl_tol,
ft_ref=ft_ref)
log.info("requires dt < {:.3e}".format(dt_tol))
@ -468,7 +484,7 @@ class StocProc_TanhSinh(_abcStocProc):
N = int(np.ceil(t_max/dt_tol))+1
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)
wmax = method_ft.find_integral_boundary(spectral_density, tol=intgr_tol/4, ref_val=1, max_val=1e6, x0=0.777)
log.info("wmax:{}".format(wmax))
tau = np.linspace(0, t_max, 35)
@ -476,7 +492,7 @@ class StocProc_TanhSinh(_abcStocProc):
d = intgr_tol+1
while d > intgr_tol:
n *= 2
I = method_fft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
I = method_ft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
x_max=wmax,
n=n,
tau_l=tau)
@ -487,7 +503,7 @@ class StocProc_TanhSinh(_abcStocProc):
tau = np.linspace(0, (N-1)*dt_tol, N)
log.info("perform numeric check of entire time axis [{},{}] n:{}".format(0, (N-1)*dt_tol, N))
num_FT = method_fft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
num_FT = method_ft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
x_max=wmax,
n=n,
tau_l=tau)
@ -511,7 +527,7 @@ class StocProc_TanhSinh(_abcStocProc):
plt.show()
assert d <= intgr_tol, "d:{}, intgr_tol:{}".format(d, intgr_tol)
yk, wk = method_fft.get_x_w_and_dt(n, wmax)
yk, wk = method_ft.get_x_w_and_dt(n, wmax)
self.omega_k = yk
self.fl = np.sqrt(wk*spectral_density(self.omega_k)/np.pi)