mirror of
https://github.com/vale981/stocproc
synced 2025-03-05 09:41:42 -05:00
testing, cleaning and docs
This commit is contained in:
Richard Hartmann
parent
9f6bdbbc62
commit
8c117e4468
7 changed files with 398 additions and 420 deletions
|
|
@ -1,4 +1,4 @@
|
|||
StocProc_FFT_tol
|
||||
----------------
|
||||
StocProc_FFT
|
||||
------------
|
||||
|
||||
.. autoclass:: stocproc.stocproc.StocProc_FFT_tol
|
||||
.. autoclass:: stocproc.stocproc.StocProc_FFT
|
||||
|
|
@ -1,9 +1,18 @@
|
|||
StocProc_KLE_tol
|
||||
----------------
|
||||
StocProc_KLE
|
||||
------------
|
||||
|
||||
|
||||
.. py:currentmodule:: stocproc
|
||||
|
||||
.. autoclass:: StocProc_KLE_tol
|
||||
.. autoclass:: StocProc_KLE
|
||||
:members:
|
||||
:inherited-members:
|
||||
:undoc-members:
|
||||
|
||||
.. 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
|
||||
|
|
|
|||
|
|
@ -56,6 +56,5 @@ import sys
|
|||
if sys.version_info.major < 3:
|
||||
raise SystemError("no support for Python 2")
|
||||
|
||||
from .stocproc import StocProc_FFT_tol
|
||||
from .stocproc import StocProc_FFT
|
||||
from .stocproc import StocProc_KLE
|
||||
from .stocproc import StocProc_KLE_tol
|
||||
|
|
@ -10,13 +10,13 @@ from . import tools
|
|||
import logging
|
||||
log = logging.getLogger(__name__)
|
||||
|
||||
def solve_hom_fredholm(r, w, eig_val_min=None):
|
||||
def solve_hom_fredholm(r, w):
|
||||
r"""Solves the discrete homogeneous Fredholm equation of the second kind
|
||||
|
||||
.. math:: \int_0^{t_\mathrm{max}} \mathrm{d}s R(t-s) u(s) = \lambda u(t)
|
||||
|
||||
Quadrature approximation of the integral gives a discrete representation of the
|
||||
basis functions and leads to a regular eigenvalue problem.
|
||||
Quadrature approximation of the integral gives a discrete representation
|
||||
which leads to a regular eigenvalue problem.
|
||||
|
||||
.. math:: \sum_i w_i R(t_j-s_i) u(s_i) = \lambda u(t_j) \equiv \mathrm{diag(w_i)} \cdot R \cdot u = \lambda u
|
||||
|
||||
|
|
@ -27,18 +27,31 @@ def solve_hom_fredholm(r, w, eig_val_min=None):
|
|||
|
||||
.. math:: D \cdot R \cdot D \cdot D \cdot u = \lambda D \cdot u \equiv \tilde R \tilde u = \lambda \tilde u
|
||||
|
||||
where :math:`\tilde R` is hermitian and :math:`u = D^{-1}\tilde u`
|
||||
|
||||
Setting :math:`t_i = s_i` and due to the definition of the correlation function the
|
||||
matrix :math:`r_{ij} = R(t_i, s_j)` is hermitian.
|
||||
where :math:`\tilde R` is hermitian and :math:`u = D^{-1}\tilde u`.
|
||||
|
||||
:param r: correlation matrix :math:`R(t_j-s_i)`
|
||||
:param w: integrations weights :math:`w_i`
|
||||
(they have to correspond to the discrete time :math:`t_i`)
|
||||
:param eig_val_min: discards all eigenvalues and eigenvectos with
|
||||
:math:`\lambda_i < \mathtt{eig\_val\_min} / \mathrm{max}(\lambda)`
|
||||
|
||||
|
||||
:return: eigenvalues, eigenvectos (eigenvectos are stored in the normal numpy fashion, ordered in decreasing order)
|
||||
|
||||
.. seealso::
|
||||
There are various convenient functions to calculate the integration weights and times
|
||||
to approximate the integral over the interval [0, t_max] using ng points.
|
||||
|
||||
:py:func:`get_mid_point_weights_times`,
|
||||
:py:func:`get_trapezoidal_weights_times`,
|
||||
:py:func:`get_simpson_weights_times`,
|
||||
:py:func:`get_four_point_weights_times`,
|
||||
:py:func:`get_gauss_legendre_weights_times`,
|
||||
:py:func:`get_tanh_sinh_weights_times`
|
||||
|
||||
.. note::
|
||||
It has been noticed that the performance of the various weights depends on the auto correlation
|
||||
function. As default one should use the simpson weights. 'four point', 'gauss legendre' and 'tanh sinh'
|
||||
might perform better for auto correlation function that decay slowly. Their advantage becomes evident
|
||||
for a large numbers of grid points only. So if one cares about relative differences below 1e-4
|
||||
the more sophisticated weights are suitable.
|
||||
"""
|
||||
t0 = time.time()
|
||||
# weighted matrix r due to quadrature weights
|
||||
|
|
@ -47,16 +60,10 @@ def solve_hom_fredholm(r, w, eig_val_min=None):
|
|||
r = w_sqrt.reshape(n,1) * r * w_sqrt.reshape(1,n)
|
||||
eig_val, eig_vec = scipy_eigh(r, overwrite_a=True) # eig_vals in ascending
|
||||
|
||||
if eig_val_min is None:
|
||||
min_idx = 0
|
||||
else:
|
||||
min_idx = sum(eig_val < eig_val_min)
|
||||
|
||||
eig_val = eig_val[min_idx:][::-1]
|
||||
eig_vec = eig_vec[:, min_idx:][:, ::-1]
|
||||
log.debug("use {} / {} eigenfunctions".format(len(eig_val), len(w)))
|
||||
eig_val = eig_val[::-1]
|
||||
eig_vec = eig_vec[:, ::-1]
|
||||
eig_vec = np.reshape(1/w_sqrt, (n,1)) * eig_vec
|
||||
log.debug("discrete fredholm equation of size {} solved [{:.2e}]".format(n, time.time() - t0))
|
||||
log.debug("discrete fredholm equation of size {} solved [{:.2e}s]".format(n, time.time() - t0))
|
||||
return eig_val, eig_vec
|
||||
|
||||
def align_eig_vec(eig_vec):
|
||||
|
|
@ -89,72 +96,72 @@ def _calc_corr_matrix(s, bcf):
|
|||
r[:, i] = bcf_n[idx:idx + n_]
|
||||
return r
|
||||
|
||||
def opt_fredholm(kernel, lam, t, u, meth, ntimes):
|
||||
def fopt(x, p, k):
|
||||
n = len(x)//2 + 1
|
||||
lam = x[0]
|
||||
u = x[1:n] + 1j*x[n:]
|
||||
return np.log10((np.sum(np.abs(k.dot(u) - lam * u) ** p)) ** (1 / p))
|
||||
|
||||
for i in range(ntimes):
|
||||
print("iter {}/{} lam {}".format(i+1, ntimes, lam))
|
||||
u_intp = tools.ComplexInterpolatedUnivariateSpline(t, u)
|
||||
ng = 2*(len(t)-1)+1
|
||||
t, w = meth(t[-1], ng)
|
||||
u_new = u_intp(t)
|
||||
k_new = kernel(t.reshape(-1,1)-t.reshape(1,-1))*w.reshape(-1,1)
|
||||
p = 5
|
||||
r = minimize(fopt, x0=np.hstack((lam, u_new.real, u_new.imag)), args=(p, k_new), method='CG')
|
||||
x = r.x
|
||||
n = len(x) // 2 + 1
|
||||
lam = x[0]
|
||||
u = x[1:n] + 1j * x[n:]
|
||||
print("func val {} ng {}".format(r.fun, ng))
|
||||
|
||||
return lam, u, t
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def get_mid_point_weights_times(t_max, num_grid_points):
|
||||
r"""Returns the time gridpoints and wiehgts for numeric integration via **mid point rule**.
|
||||
|
||||
The idea is to use :math:`N_\mathrm{GP}` time grid points located in the middle
|
||||
each of the :math:`N_\mathrm{GP}` equally distributed subintervals of :math:`[0, t_\mathrm{max}]`.
|
||||
|
||||
.. math:: t_i = \left(i + \frac{1}{2}\right)\frac{t_\mathrm{max}}{N_\mathrm{GP}} \qquad i = 0,1, ... N_\mathrm{GP} - 1
|
||||
r"""Returns the weights and grid points for numeric integration via **mid point rule**.
|
||||
|
||||
The corresponding trivial weights for integration are
|
||||
|
||||
.. math:: w_i = \Delta t = \frac{t_\mathrm{max}}{N_\mathrm{GP}} \qquad i = 0,1, ... N_\mathrm{GP} - 1
|
||||
|
||||
:param t_max: end of the interval for the time grid :math:`[0,t_\mathrm{max}]`
|
||||
(note: this would corespond to an integration from :math:`0-\Delta t / 2`
|
||||
to :math:`t_\mathrm{max}+\Delta t /2`)
|
||||
:param num_grid_points: number of
|
||||
:param num_grid_points: number of grid points N
|
||||
:return: location of the grid points, corresponding weights
|
||||
|
||||
Because this function is intended to provide the weights to be used in :py:func:`solve_hom_fredholm`
|
||||
is stretches the homogeneous weights over the grid points starting from 0 up to t_max, so the
|
||||
term min_point is somewhat miss leading. This is possible because we want to simulate
|
||||
stationary stochastic processes which allows :math:`\alpha(t_i+\Delta - (s_j+\Delta)) = \alpha(t_i-s_j)`.
|
||||
|
||||
The N grid points are located at
|
||||
|
||||
.. math:: t_i = i \frac{t_\mathrm{max}}{N-1} \qquad i = 0,1, ... N - 1
|
||||
|
||||
and the corresponding weights are
|
||||
|
||||
.. math:: w_i = \Delta t = \frac{t_\mathrm{max}}{N-1} \qquad i = 0,1, ... N - 1
|
||||
"""
|
||||
# generate mid points
|
||||
t, delta_t = np.linspace(0, t_max, num_grid_points, retstep = True)
|
||||
# equal weights for grid points
|
||||
w = np.ones(num_grid_points)*delta_t
|
||||
return t, w
|
||||
|
||||
def get_trapezoidal_weights_times(t_max, num_grid_points):
|
||||
# generate mid points
|
||||
r"""Returns the weights and grid points for numeric integration via **trapezoidal rule**.
|
||||
|
||||
:param t_max: end of the interval for the time grid :math:`[0,t_\mathrm{max}]`
|
||||
:param num_grid_points: number of grid points N
|
||||
:return: location of the grid points, corresponding weights
|
||||
|
||||
The N grid points are located at
|
||||
|
||||
.. math:: t_i = i \frac{t_\mathrm{max}}{N-1} \qquad i = 0,1, ... N - 1
|
||||
|
||||
and the corresponding weights are
|
||||
|
||||
.. math:: w_0 = w_{N-1} = \Delta t /2 \qquad w_i = \Delta t \quad 0 < i < N - 1
|
||||
\qquad \Delta t = \frac{t_\mathrm{max}}{N-1}
|
||||
"""
|
||||
t, delta_t = np.linspace(0, t_max, num_grid_points, retstep = True)
|
||||
# equal weights for grid points
|
||||
w = np.ones(num_grid_points)*delta_t
|
||||
w[0] /= 2
|
||||
w[-1] /= 2
|
||||
return t, w
|
||||
|
||||
def get_simpson_weights_times(t_max, num_grid_points):
|
||||
r"""Returns the weights and grid points for numeric integration via **simpson rule**.
|
||||
|
||||
:param t_max: end of the interval for the time grid :math:`[0,t_\mathrm{max}]`
|
||||
:param num_grid_points: number of grid points N (needs to be odd)
|
||||
:return: location of the grid points, corresponding weights
|
||||
|
||||
The N grid points are located at
|
||||
|
||||
.. math:: t_i = i \frac{t_\mathrm{max}}{N-1} \qquad i = 0,1, ... N - 1
|
||||
|
||||
and the corresponding weights are
|
||||
|
||||
.. math:: w_0 = w_{N-1} = 1/3 \Delta t \qquad w_{2i} = 2/3 \Delta t \quad 0 < i < (N - 1)/2
|
||||
.. math:: \qquad w_{2i+1} = 4/3 \Delta t \quad 0 \leq i < (N - 1)/2 \qquad \Delta t = \frac{t_\mathrm{max}}{N-1}
|
||||
"""
|
||||
if num_grid_points % 2 != 1:
|
||||
raise RuntimeError("simpson weights needs grid points ng such that ng = 2*k+1, but git ng={}".format(num_grid_points))
|
||||
# generate mid points
|
||||
t, delta_t = np.linspace(0, t_max, num_grid_points, retstep = True)
|
||||
# equal weights for grid points
|
||||
w = np.empty(num_grid_points, dtype=np.float64)
|
||||
w[0::2] = 2/3*delta_t
|
||||
w[1::2] = 4/3*delta_t
|
||||
|
|
@ -163,30 +170,76 @@ def get_simpson_weights_times(t_max, num_grid_points):
|
|||
return t, w
|
||||
|
||||
def get_four_point_weights_times(t_max, num_grid_points):
|
||||
r"""Returns the weights and grid points for numeric integration via **four-point Newton-Cotes rule**.
|
||||
|
||||
:param t_max: end of the interval for the time grid :math:`[0,t_\mathrm{max}]`
|
||||
:param num_grid_points: number of grid points N (needs to be (4k+1) where k is an integer greater 0)
|
||||
:return: location of the grid points, corresponding weights
|
||||
|
||||
The N grid points are located at
|
||||
|
||||
.. math:: t_i = i \frac{t_\mathrm{max}}{N-1} \qquad i = 0,1, ... N - 1
|
||||
|
||||
and the corresponding weights are
|
||||
|
||||
.. math:: w_0 = w_{N-1} = 28/90 \Delta t
|
||||
.. math:: w_{4i+1} = w_{4i+3} = 128/90 \Delta t \quad 0 \leq i < (N - 1)/4
|
||||
.. math:: w_{4i+2} = 48/90 \Delta t \quad 0 \leq i < (N - 1)/4
|
||||
.. math:: w_{4i} = 56/90 \Delta t \quad 0 < i < (N - 1)/4
|
||||
.. math:: \Delta t = \frac{t_\mathrm{max}}{N-1}
|
||||
"""
|
||||
if num_grid_points % 4 != 1:
|
||||
raise RuntimeError("four point weights needs grid points ng such that ng = 4*k+1, but git ng={}".format(num_grid_points))
|
||||
# generate mid points
|
||||
t, delta_t = np.linspace(0, t_max, num_grid_points, retstep = True)
|
||||
# equal weights for grid points
|
||||
w = np.empty(num_grid_points, dtype=np.float64)
|
||||
w[0::4] = 2 * 7 * 4 / 90 * delta_t
|
||||
w[1::4] = 1 * 32 * 4 / 90 * delta_t
|
||||
w[2::4] = 1 * 12 * 4 / 90 * delta_t
|
||||
w[3::4] = 1 * 32 * 4 / 90 * delta_t
|
||||
w[0] = 1 * 7 * 4 / 90 * delta_t
|
||||
w[-1] = 1 * 7 * 4 / 90 * delta_t
|
||||
w[0::4] = 56 / 90 * delta_t
|
||||
w[1::4] = 128 / 90 * delta_t
|
||||
w[2::4] = 48 / 90 * delta_t
|
||||
w[3::4] = 128 / 90 * delta_t
|
||||
w[0] = 28 / 90 * delta_t
|
||||
w[-1] = 28 / 90 * delta_t
|
||||
return t, w
|
||||
|
||||
def get_gauss_legendre_weights_times(t_max, num_grid_points):
|
||||
"""Returns the weights and grid points for numeric integration via **Gauss integration**
|
||||
by expanding the function in terms of Legendre Polynomials.
|
||||
|
||||
:param t_max: end of the interval for the time grid :math:`[0,t_\mathrm{max}]`
|
||||
:param num_grid_points: number of grid points N
|
||||
:return: location of the grid points, corresponding weights
|
||||
"""
|
||||
return gquad.gauss_nodes_weights_legendre(n = num_grid_points, low=0, high=t_max)
|
||||
|
||||
def get_sinh_tanh_weights_times(t_max, num_grid_points):
|
||||
"""
|
||||
inspired by Tanh-Sinh High-Precision Quadrature - David H. Bailey
|
||||
def get_tanh_sinh_weights_times(t_max, num_grid_points):
|
||||
r"""Returns the weights and grid points for numeric integration via **Tanh-Sinh integration**. The idea is to
|
||||
transform the integral over a finite interval :math:`x \in [-1, 1]` via the variable transformation
|
||||
|
||||
.. math :: x = \tanh(\pi/2 \sinh(t))
|
||||
|
||||
to a integral over the entire real axis :math:`t \in [-\infty,\infty]` but where the new
|
||||
transformed integrand decay rapidly such that a simply midpoint rule performs very well.
|
||||
|
||||
inspired by 'Tanh-Sinh High-Precision Quadrature - David H. Bailey'
|
||||
|
||||
:param t_max: end of the interval for the time grid :math:`[0,t_\mathrm{max}]`
|
||||
:param num_grid_points: number of grid points N (needs to be odd)
|
||||
:return: location of the grid points, corresponding weights
|
||||
|
||||
For a fixed small parameter h the location of the grid points read
|
||||
|
||||
.. math :: x_i = \tanh(\pi/2 \sinh(ih)
|
||||
|
||||
with corresponding weights
|
||||
|
||||
.. math :: w_i = \frac{\pi/2 \cosh(ih)}{\cosh^2(\pi/2 \sinh(ih))}
|
||||
|
||||
where i can be any integer. For a given number of grid points N, h is chosen such that
|
||||
:math:`w_{(N-1)/2} < 10^{-14}` which implies odd N. With that particular h :math:`x_i` and
|
||||
:math:`w_i` are calculated for :math:`-(N-1)/2 \leq i \leq (N-1)/2`. Afterwards the :math:`x_i` are linearly
|
||||
scaled such that :math:`x_{-(N-1)/2} = 0` and :math:`x_{(N-1)/2} = t_\mathrm{max}`.
|
||||
"""
|
||||
def get_h_of_N(N):
|
||||
"""returns the stepsize h for sinh tanh quad for a given number of points N
|
||||
r"""returns the stepsize h for sinh tanh quad for a given number of points N
|
||||
such that the smallest weights are about 1e-14
|
||||
"""
|
||||
a = 16.12087683080651
|
||||
|
|
@ -210,18 +263,6 @@ def get_sinh_tanh_weights_times(t_max, num_grid_points):
|
|||
t = np.hstack((y_plus[-1:0:-1], (2-y_plus)))*t_max/2
|
||||
return t, w
|
||||
|
||||
|
||||
_WC_ = 2
|
||||
_S_ = 0.6
|
||||
_GAMMA_S_PLUS_1 = gamma(_S_ + 1)
|
||||
def lac(t):
|
||||
"""lorenzian bath correlation function"""
|
||||
return np.exp(- np.abs(t) - 1j * _WC_ * t)
|
||||
|
||||
def oac(t):
|
||||
"""ohmic bath correlation function"""
|
||||
return (1 + 1j*(t))**(-(_S_+1)) * _GAMMA_S_PLUS_1 / np.pi
|
||||
|
||||
def get_rel_diff(corr, t_max, ng, weights_times, ng_fac):
|
||||
t, w = weights_times(t_max, ng)
|
||||
r = corr(t.reshape(-1, 1) - t.reshape(1, -1))
|
||||
|
|
@ -252,14 +293,69 @@ def get_rel_diff(corr, t_max, ng, weights_times, ng_fac):
|
|||
return tfine, rd
|
||||
|
||||
def auto_ng(corr, t_max, ngfac=2, meth=get_mid_point_weights_times, tol=1e-3, diff_method='full', dm_random_samples=10**4):
|
||||
r"""increase the number of gridpoints until the desired accuracy is met
|
||||
|
||||
This function increases the number of grid points of the discrete Fredholm equation exponentially until
|
||||
a given accuracy is met. The accuracy is determined from the deviation of the approximated
|
||||
auto correlation of the Karhunen-Loève expansion from the given reference auto correlation.
|
||||
|
||||
.. math::
|
||||
|
||||
\Delta(n) = \max_{t,s \in [0,t_\mathrm{max}]}\left( \Big | \alpha(t-s) - \sum_{i=1}^n \lambda_i u_i(t) u_i^\ast(s) \Big | \right )
|
||||
|
||||
The procedure works as follows:
|
||||
1) Solve the discrete Fredholm equation on a grid with ng points.
|
||||
This gives ng eigenvalues/vectors where each ng-dimensional vector approximates the continuous eigenfunction.
|
||||
(:math:`t, u_i(t) \leftrightarrow t_k, u_{ik}` where the :math:`t_k` depend on the integration weights
|
||||
method)
|
||||
2) Approximate the eigenfunction on a finer, equidistant grid
|
||||
(:math:`ng_\mathrm{fine} = ng_\mathrm{fac}(ng-1)+1`) using
|
||||
|
||||
.. math::
|
||||
|
||||
u_i(t) = \frac{1}{\lambda_i} \int_0^{t_\mathrm{max}} \mathrm{d}s \; \alpha(t-s) u_i(s)
|
||||
\approx \frac{1}{\lambda_i} \sum_k w_k \alpha(t-s_k) u_{ik}
|
||||
|
||||
According to the Numerical Recipes [1] this interpolation should perform better that simple
|
||||
spline interpolation. However it turns that this is not the case in general (e.g. for exponential
|
||||
auto correlation functions the spline interpolation performs better). For that reason it might be
|
||||
usefull to set ngfac to 1 which will skip the integral interpolation
|
||||
3) Use the eigenfunction on the fine grid to setup a cubic spline interpolation.
|
||||
4) Use the spline interpolation to estimate the deviation :math:`\Delta(n)`. When using diff_method = 'full'
|
||||
the maximization is performed over all :math:`t'_i, s'_j` where :math:`t'_i = (t_i + t_{i+1})/2` and
|
||||
:math:`s'_i = (s_i + s_{i+1})/2` with :math:`i,j = 0, \, ...\, , ng_\mathrm{fine}-2`. It is expected that
|
||||
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
|
||||
which is needed for the 'simpson' and 'fourpoint' integration weights)
|
||||
|
||||
:param corr: the auto correlation function
|
||||
:param t_max: specifies the interval [0, t_max] for which the stochastic process can be evaluated
|
||||
:param ngfac: specifies the fine grid to use for the spline interpolation, the intermediate points are
|
||||
calculated using integral interpolation
|
||||
:param meth: the method for calculation integration weights and times, a callable or one of the following strings
|
||||
'midpoint' ('midp'), 'trapezoidal' ('trapz'), 'simpson' ('simp'), 'fourpoint' ('fp'),
|
||||
'gauss_legendre' ('gl'), 'tanh_sinh' ('ts')
|
||||
:param tol: defines the success criterion max(abs(corr_exact - corr_reconstr)) < tol
|
||||
:param diff_method: either 'full' or 'random', determines the points where the above success criterion is evaluated,
|
||||
'full': full grid in between the fine grid, such that the spline interpolation error is expected to be maximal
|
||||
'random': pick a fixed number of random times t and s within the interval [0, t_max]
|
||||
:param dm_random_samples: the number of random times used for diff_method 'random'
|
||||
:return: an array containing the necessary eigenfunctions of the Karhunen-Loève expansion for sampling the
|
||||
stochastic processes
|
||||
|
||||
[1] 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. (pp. 990)
|
||||
"""
|
||||
time_start = time.time()
|
||||
if diff_method == 'full':
|
||||
pass
|
||||
elif diff_method == 'random':
|
||||
t_rand = np.random.rand(dm_random_samples) * t_max
|
||||
s_rand = np.random.rand(dm_random_samples) * t_max
|
||||
alpha_ref = corr(t_rand - s_rand)
|
||||
diff = -alpha_ref
|
||||
else:
|
||||
raise ValueError("unknown diff_method '{}', use 'full' or 'random'".format(diff_method))
|
||||
alpha_0 = np.abs(corr(0))
|
||||
|
|
@ -271,77 +367,107 @@ def auto_ng(corr, t_max, ngfac=2, meth=get_mid_point_weights_times, tol=1e-3, di
|
|||
time_spline = 0
|
||||
time_calc_diff = 0
|
||||
|
||||
if isinstance(meth, str):
|
||||
meth = str_meth_to_meth(meth)
|
||||
|
||||
k = 4
|
||||
while True:
|
||||
k += 1
|
||||
ng = 2 ** k + 1
|
||||
log.debug("check {} grid points".format(ng))
|
||||
log.info("check {} grid points".format(ng))
|
||||
t, w = meth(t_max, ng)
|
||||
t0 = time.time()
|
||||
r = _calc_corr_matrix(t, corr)
|
||||
|
||||
t0 = time.time() # efficient way to construct the
|
||||
r = _calc_corr_matrix(t, corr) # auto correlation matrix r
|
||||
time_calc_ac += (time.time() - t0)
|
||||
|
||||
t0 = time.time()
|
||||
_eig_val, _eig_vec = solve_hom_fredholm(r, w)
|
||||
t0 = time.time() # solve the dicrete fredholm equation
|
||||
_eig_val, _eig_vec = solve_hom_fredholm(r, w) # using integration weights w
|
||||
time_fredholm += (time.time() - t0)
|
||||
|
||||
tfine, dt = np.linspace(0, t_max, (ng - 1) * ngfac + 1, retstep=True)
|
||||
tsfine = tfine[:-1] + dt/2
|
||||
tfine, dt = np.linspace(0, t_max, (ng - 1) * ngfac + 1, retstep=True) # setup fine
|
||||
tsfine = tfine[:-1] + dt/2 # and super fine time grid
|
||||
|
||||
t0 = time.time()
|
||||
alpha_k = _calc_corr_min_t_plus_t(tfine, corr)
|
||||
time_calc_ac += (time.time() - t0)
|
||||
t0 = time.time() # efficient way to calculate the auto correlation
|
||||
alpha_k = _calc_corr_min_t_plus_t(tfine, corr) # from -tmax untill tmax on the fine grid
|
||||
time_calc_ac += (time.time() - t0) # needed for integral interpolation
|
||||
|
||||
if diff_method == 'full':
|
||||
ng_sfine = len(tsfine)
|
||||
alpha_ref = np.empty(shape=(ng_sfine, ng_sfine), dtype=np.complex128)
|
||||
for i in range(ng_sfine):
|
||||
idx = ng_sfine - 1 - i
|
||||
idx = ng_sfine - i
|
||||
alpha_ref[:, i] = alpha_k[idx:idx + ng_sfine] # note we can use alpha_k as
|
||||
# alpha(ti+dt/2 - (tj+dt/2)) = alpha(ti - tj)
|
||||
diff = -alpha_ref
|
||||
old_md = np.inf
|
||||
ui_fine_all = []
|
||||
sqrt_lambda_ui_fine_all = []
|
||||
for i in range(ng):
|
||||
evec = _eig_vec[:, i]
|
||||
sqrt_eval = np.sqrt(_eig_val[i])
|
||||
if ngfac != 1:
|
||||
t0 = time.time()
|
||||
ui_fine = stocproc_c.eig_func_interp(delta_t_fac=ngfac,
|
||||
time_axis=t,
|
||||
alpha_k=alpha_k,
|
||||
weights=w,
|
||||
eigen_val=sqrt_eval,
|
||||
eigen_vec=evec)
|
||||
# when using sqrt_lambda instead of lambda we get sqrt_lamda time u
|
||||
# which is the quantity needed for the stochastic process
|
||||
# generation
|
||||
sqrt_lambda_ui_fine = stocproc_c.eig_func_interp(delta_t_fac=ngfac,
|
||||
time_axis=t,
|
||||
alpha_k=alpha_k,
|
||||
weights=w,
|
||||
eigen_val=sqrt_eval,
|
||||
eigen_vec=evec)
|
||||
time_integr_intp += (time.time() - t0)
|
||||
else:
|
||||
ui_fine = evec
|
||||
ui_fine_all.append(ui_fine)
|
||||
sqrt_lambda_ui_fine = evec*sqrt_eval
|
||||
|
||||
sqrt_lambda_ui_fine_all.append(sqrt_lambda_ui_fine)
|
||||
|
||||
# setup cubic spline interpolator
|
||||
t0 = time.time()
|
||||
ui_spl = tools.ComplexInterpolatedUnivariateSpline(tfine, ui_fine)
|
||||
sqrt_lambda_ui_spl = tools.ComplexInterpolatedUnivariateSpline(tfine, sqrt_lambda_ui_fine)
|
||||
time_spline += (time.time() - t0)
|
||||
|
||||
# calculate the max deviation
|
||||
t0 = time.time()
|
||||
if diff_method == 'random':
|
||||
ui_t = ui_spl(t_rand)
|
||||
ui_s = ui_spl(s_rand)
|
||||
ui_t = sqrt_lambda_ui_spl(t_rand)
|
||||
ui_s = sqrt_lambda_ui_spl(s_rand)
|
||||
diff += ui_t * np.conj(ui_s)
|
||||
elif diff_method == 'full':
|
||||
ui_super_fine = ui_spl(tsfine)
|
||||
ui_super_fine = sqrt_lambda_ui_spl(tsfine)
|
||||
diff += ui_super_fine.reshape(-1, 1) * np.conj(ui_super_fine.reshape(1, -1))
|
||||
md = np.max(np.abs(diff)) / alpha_0
|
||||
time_calc_diff += (time.time() - t0)
|
||||
|
||||
log.debug("num evec {} -> max diff {:.3e}".format(i+1, md))
|
||||
if old_md < md:
|
||||
log.debug("max diff increased -> break, use higher ng")
|
||||
break
|
||||
old_md = md
|
||||
#if old_md < md:
|
||||
# log.info("max diff increased -> break, use higher ng")
|
||||
# break
|
||||
#old_md = md
|
||||
if md < tol:
|
||||
time_total = time_calc_diff + time_spline + time_integr_intp + time_calc_ac + time_fredholm
|
||||
log.debug("calc_ac {:.3%}, fredholm {:.3%}, integr_intp {:.3%}, spline {:.3%}, calc_diff {:.3%}".format(
|
||||
time_calc_ac/time_total, time_fredholm/time_total, time_integr_intp/time_total,
|
||||
time_spline/time_total, time_calc_diff/time_total))
|
||||
log.debug("auto ng SUCCESSFUL max diff {:.3e} < tol {:.3e} ng {} num evec {}".format(md, tol, ng, i+1))
|
||||
return np.asarray(ui_fine_all)
|
||||
time_overall = time.time() - time_start
|
||||
time_rest = time_overall - time_total
|
||||
|
||||
log.info("calc_ac {:.3%}, fredholm {:.3%}, integr_intp {:.3%}, spline {:.3%}, calc_diff {:.3%}, rest {:.3%}".format(
|
||||
time_calc_ac/time_overall, time_fredholm/time_overall, time_integr_intp/time_overall,
|
||||
time_spline/time_overall, time_calc_diff/time_overall, time_rest/time_overall))
|
||||
log.info("auto ng SUCCESSFUL max diff {:.3e} < tol {:.3e} ng {} num evec {}".format(md, tol, ng, i+1))
|
||||
return np.asarray(sqrt_lambda_ui_fine_all)
|
||||
log.info("ng {} yields md {:.3e}".format(ng, md))
|
||||
|
||||
def str_meth_to_meth(meth):
|
||||
if (meth == 'midpoint') or (meth == 'midp'):
|
||||
return get_mid_point_weights_times
|
||||
elif (meth == 'trapezoidal') or (meth == 'trapz'):
|
||||
return get_trapezoidal_weights_times
|
||||
elif (meth == 'simpson') or (meth == 'simp'):
|
||||
return get_simpson_weights_times
|
||||
elif (meth == 'fourpoint') or (meth == 'fp'):
|
||||
return get_four_point_weights_times
|
||||
elif (meth == 'gauss_legendre') or (meth == 'gl'):
|
||||
return get_gauss_legendre_weights_times
|
||||
elif (meth == 'tanh_sinh') or (meth == 'ts'):
|
||||
return get_tanh_sinh_weights_times
|
||||
else:
|
||||
raise ValueError("unknown method to get integration weights '{}'".format(meth))
|
||||
|
|
@ -24,6 +24,7 @@ associated spectral density. The number scales linear with the time interval of
|
|||
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.
|
||||
"""
|
||||
import abc
|
||||
import numpy as np
|
||||
import time
|
||||
|
||||
|
|
@ -35,7 +36,7 @@ from .tools import ComplexInterpolatedUnivariateSpline
|
|||
import logging
|
||||
log = logging.getLogger(__name__)
|
||||
|
||||
class _absStocProc(object):
|
||||
class _absStocProc(abc.ABC):
|
||||
r"""
|
||||
Abstract base class to stochastic process interface
|
||||
|
||||
|
|
@ -49,20 +50,18 @@ class _absStocProc(object):
|
|||
|
||||
|
||||
"""
|
||||
def __init__(self, t_max, num_grid_points, seed=None, k=3):
|
||||
def __init__(self, t_max, num_grid_points, seed=None):
|
||||
r"""
|
||||
:param t_max: specify time axis as [0, t_max]
|
||||
:param num_grid_points: number of equidistant times on that axis
|
||||
:param seed: if not ``None`` set seed to ``seed``
|
||||
:param verbose: 0: no output, 1: informational output, 2: (eventually some) debug info
|
||||
:param k: polynomial degree used for spline interpolation
|
||||
"""
|
||||
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
|
||||
self._k = k
|
||||
self._seed = seed
|
||||
if seed is not None:
|
||||
np.random.seed(seed)
|
||||
|
|
@ -83,10 +82,11 @@ class _absStocProc(object):
|
|||
else:
|
||||
if self._interpolator is None:
|
||||
t0 = time.time()
|
||||
self._interpolator = ComplexInterpolatedUnivariateSpline(self.t, self._z, k=self._k)
|
||||
self._interpolator = ComplexInterpolatedUnivariateSpline(self.t, self._z, k=3)
|
||||
log.debug("created interpolator [{:.2e}s]".format(time.time() - t0))
|
||||
return self._interpolator(t)
|
||||
|
||||
|
||||
@abc.abstractmethod
|
||||
def _calc_z(self, y):
|
||||
r"""
|
||||
maps the normal distributed complex valued random variables y to the stochastic process
|
||||
|
|
@ -94,7 +94,8 @@ class _absStocProc(object):
|
|||
:return: the stochastic process, array of complex numbers
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
@abc.abstractmethod
|
||||
def get_num_y(self):
|
||||
r"""
|
||||
:return: number of complex random variables needed to calculate the stochastic process
|
||||
|
|
@ -136,127 +137,8 @@ class _absStocProc(object):
|
|||
self._z = self._calc_z(y)
|
||||
log.debug("proc_cnt:{} new process generated [{:.2e}s]".format(self._proc_cnt, time.time() - t0))
|
||||
|
||||
METHOD = 'midp'
|
||||
|
||||
class StocProc_KLE(_absStocProc):
|
||||
r"""
|
||||
class to simulate stochastic processes using KLE method
|
||||
- Solve fredholm equation on grid with ``ng_fredholm nodes`` (trapezoidal_weights).
|
||||
If case ``ng_fredholm`` is ``None`` set ``ng_fredholm = num_grid_points``. In general it should
|
||||
hold ``ng_fredholm < num_grid_points`` and ``num_grid_points = 10*ng_fredholm`` might be a good ratio.
|
||||
- Calculate discrete stochastic process (using interpolation solution of fredholm equation) with num_grid_points nodes
|
||||
- invoke spline interpolator when calling
|
||||
"""
|
||||
def __init__(self, r_tau, t_max, ng_fredholm, ng_fac=4, seed=None, sig_min=1e-5, k=3, align_eig_vec=False):
|
||||
r"""
|
||||
:param r_tau: auto correlation function of the stochastic process
|
||||
:param t_max: specify time axis as [0, t_max]
|
||||
:param seed: if not ``None`` set seed to ``seed``
|
||||
:param sig_min: eigenvalue threshold (see KLE method to generate stochastic processes)
|
||||
:param verbose: 0: no output, 1: informational output, 2: (eventually some) debug info
|
||||
:param k: polynomial degree used for spline interpolation
|
||||
"""
|
||||
|
||||
# this lengthy part will be skipped when init class from dump, as _A and alpha_k will be stored
|
||||
t0 = time.time()
|
||||
if METHOD == 'midp':
|
||||
t, w = method_kle.get_mid_point_weights_times(t_max, ng_fredholm)
|
||||
elif METHOD == 'simp':
|
||||
t, w = method_kle.get_simpson_weights_times(t_max, ng_fredholm)
|
||||
|
||||
r = self._calc_corr_matrix(t, r_tau)
|
||||
_eig_val, _eig_vec = method_kle.solve_hom_fredholm(r, w, sig_min ** 2)
|
||||
if align_eig_vec:
|
||||
for i in range(_eig_vec.shape[1]):
|
||||
s = np.sum(_eig_vec[:, i])
|
||||
phase = np.exp(1j * np.arctan2(np.real(s), np.imag(s)))
|
||||
_eig_vec[:, i] /= phase
|
||||
_sqrt_eig_val = np.sqrt(_eig_val)
|
||||
_A = w.reshape(-1, 1) * _eig_vec / _sqrt_eig_val.reshape(1, -1)
|
||||
ng_fine = ng_fac * (ng_fredholm - 1) + 1
|
||||
alpha_k = self._calc_corr_min_t_plus_t(s=np.linspace(0, t_max, ng_fine), bcf=r_tau)
|
||||
log.debug("new KLE StocProc class prepared [{:.2e}]".format(time.time() - t0))
|
||||
|
||||
data = (_A, alpha_k, seed, k, t_max, ng_fac)
|
||||
self.__setstate__(data)
|
||||
|
||||
# needed for getkey / getstate
|
||||
self.key = (r_tau, t_max, ng_fredholm, ng_fac, sig_min, align_eig_vec)
|
||||
|
||||
# save these guys as they are needed to estimate the autocorrelation
|
||||
self._s = t
|
||||
self._w = w
|
||||
self._eig_val = _eig_val
|
||||
self._eig_vec = _eig_vec
|
||||
|
||||
def _getkey(self):
|
||||
return self.__class__.__name__, self.key
|
||||
|
||||
def __getstate__(self):
|
||||
return self._A, self.alpha_k, self._seed, self._k, self.t_max, self.ng_fac
|
||||
|
||||
def __setstate__(self, state):
|
||||
self._A, self.alpha_k, seed, k, t_max, self.ng_fac = state
|
||||
if self.ng_fac == 1:
|
||||
self.kle_interp = False
|
||||
else:
|
||||
self.kle_interp = True
|
||||
self._one_over_sqrt_2 = 1 / np.sqrt(2)
|
||||
num_gp, self.num_y = self._A.shape
|
||||
ng_fine = self.ng_fac * (num_gp - 1) + 1
|
||||
super().__init__(t_max=t_max, num_grid_points=ng_fine, seed=seed, k=k)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@staticmethod
|
||||
def _calc_corr_min_t_plus_t(s, bcf):
|
||||
bcf_n_plus = bcf(s-s[0])
|
||||
# [bcf(-3) , bcf(-2) , bcf(-1) , bcf(0), bcf(1), bcf(2), bcf(3)]
|
||||
# == [bcf(3)^\ast, bcf(2)^\ast, bcf(1)^\ast, bcf(0), bcf(1), bcf(2), bcf(3)]
|
||||
return np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus))
|
||||
|
||||
@staticmethod
|
||||
def _calc_corr_matrix(s, bcf):
|
||||
"""calculates the matrix alpha_ij = bcf(t_i-s_j)
|
||||
|
||||
calls bcf only for s-s_0 and reconstructs the rest
|
||||
"""
|
||||
n_ = len(s)
|
||||
bcf_n = StocProc_KLE._calc_corr_min_t_plus_t(s, bcf)
|
||||
# we want
|
||||
# r = bcf(0) bcf(-1), bcf(-2)
|
||||
# bcf(1) bcf( 0), bcf(-1)
|
||||
# bcf(2) bcf( 1), bcf( 0)
|
||||
r = np.empty(shape=(n_,n_), dtype = np.complex128)
|
||||
for i in range(n_):
|
||||
idx = n_-1-i
|
||||
r[:,i] = bcf_n[idx:idx+n_]
|
||||
return r
|
||||
|
||||
def __calc_missing(self):
|
||||
raise NotImplementedError
|
||||
|
||||
def _calc_z(self, y):
|
||||
if self.kle_interp:
|
||||
_a_tmp = np.tensordot(y, self._A, axes=([0], [1]))
|
||||
_num_gp = self._A.shape[0]
|
||||
return stocproc_c.z_t(delta_t_fac = self.ng_fac,
|
||||
N1 = _num_gp,
|
||||
alpha_k = self.alpha_k,
|
||||
a_tmp = _a_tmp,
|
||||
kahanSum = True)
|
||||
|
||||
else:
|
||||
return np.tensordot(y*self._sqrt_eig_val, self._eig_vec, axes=([0], [1])).flatten()
|
||||
|
||||
def get_num_y(self):
|
||||
return self.num_y
|
||||
|
||||
|
||||
|
||||
class StocProc_KLE_tol(StocProc_KLE):
|
||||
r"""
|
||||
A class to simulate stochastic processes using Karhunen-Loève expansion (KLE) method.
|
||||
The idea is that any stochastic process can be expressed in terms of the KLE
|
||||
|
|
@ -272,7 +154,8 @@ class StocProc_KLE_tol(StocProc_KLE):
|
|||
for a given positive integral kernel :math:`R(\tau)`. It turns out that the auto correlation
|
||||
:math:`\langle Z(t)Z^\ast(s) \rangle = R(t-s)` is given by that kernel.
|
||||
|
||||
For the numeric implementation the integral equation has to be discretized
|
||||
For the numeric implementation the integral equation will be discretized
|
||||
(see :py:func:`stocproc.method_kle.solve_hom_fredholm` for details).
|
||||
|
||||
|
||||
- Solve fredholm equation on grid with ``ng_fredholm nodes`` (trapezoidal_weights).
|
||||
|
|
@ -287,7 +170,8 @@ class StocProc_KLE_tol(StocProc_KLE):
|
|||
|
||||
"""
|
||||
|
||||
def __init__(self, r_tau, t_max, tol=1e-2, ng_fac=4, seed=None, k=3, align_eig_vec=False):
|
||||
def __init__(self, r_tau, t_max, ng_fac=4, meth='fourpoint', tol=1e-2, diff_method='full', dm_random_samples=10**4,
|
||||
seed=None, align_eig_vec=False):
|
||||
"""this is init
|
||||
|
||||
:param r_tau:
|
||||
|
|
@ -298,69 +182,48 @@ class StocProc_KLE_tol(StocProc_KLE):
|
|||
:param k:
|
||||
:param align_eig_vec:
|
||||
"""
|
||||
self.tol = tol
|
||||
kwargs = {'r_tau': r_tau, 't_max': t_max, 'ng_fac': ng_fac, 'seed': seed,
|
||||
'sig_min': tol**2, 'k': k, 'align_eig_vec': align_eig_vec}
|
||||
self._auto_grid_points(**kwargs)
|
||||
# overwrite ng_fac in key from StocProc_KLE with value of tol
|
||||
# self.key = (r_tau, t_max, ng_fredholm, ng_fac, sig_min, align_eig_vec)
|
||||
self.key = (self.key[0], self.key[1], tol, self.key[3],self.key[4], self.key[5])
|
||||
|
||||
def _init_StocProc_KLE_and_get_error(self, ng, **kwargs):
|
||||
super().__init__(ng_fredholm=ng, **kwargs)
|
||||
|
||||
ng_fine = self.ng_fac*(ng-1)+1
|
||||
u_i_all_t = stocproc_c.eig_func_all_interp(delta_t_fac = self.ng_fac,
|
||||
time_axis = self._s,
|
||||
alpha_k = self.alpha_k,
|
||||
weights = self._w,
|
||||
eigen_val = self._eig_val,
|
||||
eigen_vec = self._eig_vec)
|
||||
sqrt_lambda_ui_fine = method_kle.auto_ng(corr=r_tau,
|
||||
t_max=t_max,
|
||||
ngfac=ng_fac,
|
||||
meth=meth,
|
||||
tol=tol,
|
||||
diff_method=diff_method,
|
||||
dm_random_samples=dm_random_samples)
|
||||
|
||||
u_i_all_ast_s = np.conj(u_i_all_t) #(N_gp, N_ev)
|
||||
num_ev = len(self._eig_val)
|
||||
tmp = self._eig_val.reshape(1, num_ev) * u_i_all_t #(N_gp, N_ev)
|
||||
recs_bcf = np.tensordot(tmp, u_i_all_ast_s, axes=([1],[1]))
|
||||
|
||||
refc_bcf = np.empty(shape=(ng_fine,ng_fine), dtype = np.complex128)
|
||||
for i in range(ng_fine):
|
||||
idx = ng_fine-1-i
|
||||
refc_bcf[:,i] = self.alpha_k[idx:idx+ng_fine]
|
||||
|
||||
err = np.max(np.abs(recs_bcf-refc_bcf)/np.abs(refc_bcf))
|
||||
return err
|
||||
|
||||
|
||||
def _auto_grid_points(self, **kwargs):
|
||||
err = np.inf
|
||||
c = 2
|
||||
#exponential increase to get below error threshold
|
||||
while err > self.tol:
|
||||
c *= 2
|
||||
ng = 2*c + 1
|
||||
err = self._init_StocProc_KLE_and_get_error(ng, **kwargs)
|
||||
log.info("ng {} -> err {:.3e}".format(ng, err))
|
||||
|
||||
c_low = c // 2
|
||||
c_high = c
|
||||
|
||||
while (c_high - c_low) > 1:
|
||||
c = (c_low + c_high) // 2
|
||||
ng = 2*c + 1
|
||||
err = self._init_StocProc_KLE_and_get_error(ng, **kwargs)
|
||||
log.info("ng {} -> err {:.3e}".format(ng, err))
|
||||
if err > self.tol:
|
||||
c_low = c
|
||||
else:
|
||||
c_high = c
|
||||
# inplace alignment such that re(ui(0)) >= 0 and im(ui(0)) = 0
|
||||
if align_eig_vec:
|
||||
method_kle.align_eig_vec(sqrt_lambda_ui_fine)
|
||||
|
||||
state = sqrt_lambda_ui_fine, t_max, seed
|
||||
self.__setstate__(state)
|
||||
self.key = r_tau, t_max, tol
|
||||
|
||||
def __getstate__(self):
|
||||
return self.sqrt_lambda_ui_fine, self.t_max, self._seed
|
||||
|
||||
def __setstate__(self, state):
|
||||
sqrt_lambda_ui_fine, t_max, seed = state
|
||||
num_ev, ng = sqrt_lambda_ui_fine.shape
|
||||
super().__init__(t_max = t_max,
|
||||
num_grid_points = ng,
|
||||
seed = seed)
|
||||
self.num_ev = num_ev
|
||||
self.sqrt_lambda_ui_fine = sqrt_lambda_ui_fine
|
||||
|
||||
def _calc_z(self, y):
|
||||
return np.tensordot(y, self.sqrt_lambda_ui_fine, axes=([0], [0])).flatten()
|
||||
|
||||
def get_num_y(self):
|
||||
return self.num_ev
|
||||
|
||||
|
||||
class StocProc_FFT_tol(_absStocProc):
|
||||
class StocProc_FFT(_absStocProc):
|
||||
r"""
|
||||
Simulate Stochastic Process using FFT method
|
||||
"""
|
||||
def __init__(self, spectral_density, t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2,
|
||||
seed=None, k=3, negative_frequencies=False):
|
||||
seed=None, negative_frequencies=False):
|
||||
if not negative_frequencies:
|
||||
log.info("non neg freq only")
|
||||
# assume the spectral_density is 0 for w<0
|
||||
|
|
@ -412,24 +275,22 @@ class StocProc_FFT_tol(_absStocProc):
|
|||
|
||||
super().__init__(t_max = t_max,
|
||||
num_grid_points = num_grid_points,
|
||||
seed = seed,
|
||||
k = k)
|
||||
seed = seed)
|
||||
|
||||
omega = dx*np.arange(N)
|
||||
self.yl = spectral_density(omega + 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
|
||||
self.key = bcf_ref, 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._k
|
||||
return self.yl, self.num_grid_points, self.omega_min_correction, self.t_max, self._seed
|
||||
|
||||
def __setstate__(self, state):
|
||||
self.yl, num_grid_points, self.omega_min_correction, t_max, seed, k = state
|
||||
self.yl, num_grid_points, self.omega_min_correction, t_max, seed = state
|
||||
super().__init__(t_max = t_max,
|
||||
num_grid_points = num_grid_points,
|
||||
seed = seed,
|
||||
k = k)
|
||||
|
||||
seed = seed)
|
||||
|
||||
def _calc_z(self, y):
|
||||
z = np.fft.fft(self.yl * y)[0:self.num_grid_points] * self.omega_min_correction
|
||||
|
|
|
|||
|
|
@ -482,13 +482,15 @@ def test_cython_interpolation():
|
|||
|
||||
def test_reconstr_ac():
|
||||
t_max = 15
|
||||
method_kle.auto_ng(corr=method_kle.oac,
|
||||
t_max=t_max,
|
||||
ngfac=2,
|
||||
meth=method_kle.get_mid_point_weights_times,
|
||||
tol=1e-3,
|
||||
diff_method='random',
|
||||
dm_random_samples=10**4)
|
||||
res = method_kle.auto_ng(corr=method_kle.oac,
|
||||
t_max=t_max,
|
||||
ngfac=2,
|
||||
meth=method_kle.get_mid_point_weights_times,
|
||||
tol=1e-3,
|
||||
diff_method='full',
|
||||
dm_random_samples=10 ** 4)
|
||||
print(type(res))
|
||||
print(res.shape)
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -20,11 +20,13 @@ import sys
|
|||
import os
|
||||
|
||||
import pathlib
|
||||
|
||||
p = pathlib.PosixPath(os.path.abspath(__file__))
|
||||
sys.path.insert(0, str(p.parent.parent))
|
||||
|
||||
import stocproc as sp
|
||||
import warnings
|
||||
|
||||
warnings.simplefilter('default')
|
||||
|
||||
_S_ = 0.6
|
||||
|
|
@ -33,10 +35,20 @@ _GAMMA_S_PLUS_1 = gamma(_S_ + 1)
|
|||
|
||||
def corr(tau):
|
||||
"""ohmic bath correlation function"""
|
||||
return (1 + 1j*(tau))**(-(_S_+1)) * _GAMMA_S_PLUS_1 / np.pi
|
||||
return (1 + 1j * (tau)) ** (-(_S_ + 1)) * _GAMMA_S_PLUS_1 / np.pi
|
||||
|
||||
|
||||
def spectral_density(omega):
|
||||
return omega**_S_ * np.exp(-omega)
|
||||
return omega ** _S_ * np.exp(-omega)
|
||||
|
||||
_WC_ = 2
|
||||
def lac(t):
|
||||
"""lorenzian bath correlation function"""
|
||||
return np.exp(- np.abs(t) - 1j * _WC_ * t)
|
||||
|
||||
def lsd(w):
|
||||
return 1 / (1 + (w - _WC_) ** 2)
|
||||
|
||||
|
||||
def stocproc_metatest(stp, num_samples, tol, corr, plot):
|
||||
print("generate samples")
|
||||
|
|
@ -99,26 +111,6 @@ def stocproc_metatest(stp, num_samples, tol, corr, plot):
|
|||
|
||||
|
||||
def test_stochastic_process_KLE_correlation_function(plot=False):
|
||||
"""
|
||||
generate samples using StocPrioc_KLE class
|
||||
|
||||
compare <X_t X_s> and <X_t X^ast_s> from samples with true ac functions
|
||||
|
||||
no interpolation at all
|
||||
"""
|
||||
t_max = 15
|
||||
num_grid_points = 101
|
||||
num_samples = 1000
|
||||
tol = 3e-2
|
||||
stp = sp.StocProc_KLE(r_tau = corr,
|
||||
t_max = t_max,
|
||||
ng_fredholm = num_grid_points,
|
||||
ng_fac = 4,
|
||||
seed = 0)
|
||||
stocproc_metatest(stp, num_samples, tol, corr, plot)
|
||||
|
||||
|
||||
def test_stochastic_process_KLE_tol_correlation_function(plot=False):
|
||||
"""
|
||||
generate samples using FFT method
|
||||
|
||||
|
|
@ -128,18 +120,13 @@ def test_stochastic_process_KLE_tol_correlation_function(plot=False):
|
|||
"""
|
||||
|
||||
t_max = 15
|
||||
num_samples = 1000
|
||||
num_samples = 2000
|
||||
tol = 3e-2
|
||||
stp = sp.StocProc_KLE_tol(tol = 1e-2,
|
||||
r_tau = corr,
|
||||
t_max = t_max,
|
||||
ng_fac = 4,
|
||||
seed = 0)
|
||||
stp = sp.StocProc_KLE(tol=1e-2, r_tau=corr, t_max=t_max, ng_fac=4, seed=0)
|
||||
stocproc_metatest(stp, num_samples, tol, corr, plot)
|
||||
|
||||
|
||||
|
||||
def test_stochastic_process_FFT_correlation_function(plot = False):
|
||||
def test_stochastic_process_FFT_correlation_function(plot=False):
|
||||
"""
|
||||
generate samples using FFT method
|
||||
|
||||
|
|
@ -149,26 +136,20 @@ def test_stochastic_process_FFT_correlation_function(plot = False):
|
|||
"""
|
||||
|
||||
t_max = 15
|
||||
num_samples = 1000
|
||||
num_samples = 2000
|
||||
tol = 3e-2
|
||||
stp = sp.StocProc_FFT_tol(spectral_density = spectral_density,
|
||||
t_max = t_max,
|
||||
bcf_ref = corr,
|
||||
intgr_tol = 1e-2,
|
||||
intpl_tol = 1e-2,
|
||||
seed = 0)
|
||||
stp = sp.StocProc_FFT(spectral_density=spectral_density, t_max=t_max, bcf_ref=corr, intgr_tol=1e-2, intpl_tol=1e-2,
|
||||
seed=0)
|
||||
stocproc_metatest(stp, num_samples, tol, corr, plot)
|
||||
|
||||
|
||||
def test_stocproc_dump_load():
|
||||
t_max = 15
|
||||
num_grid_points = 101
|
||||
|
||||
## STOCPROC KLE ##
|
||||
####################
|
||||
t0 = time.time()
|
||||
stp = sp.StocProc_KLE(r_tau = corr,
|
||||
t_max = t_max,
|
||||
ng_fredholm = num_grid_points,
|
||||
ng_fac = 4,
|
||||
seed = 0)
|
||||
stp = sp.StocProc_KLE(tol=1e-2, r_tau=corr, t_max=t_max, ng_fac=4, seed=0)
|
||||
t1 = time.time()
|
||||
dt1 = t1 - t0
|
||||
stp.new_process()
|
||||
|
|
@ -179,39 +160,18 @@ def test_stocproc_dump_load():
|
|||
stp2 = pickle.loads(bin_data)
|
||||
t1 = time.time()
|
||||
dt2 = t1 - t0
|
||||
assert dt2 / dt1 < 0.1 # loading should be way faster
|
||||
assert dt2 / dt1 < 0.1 # loading should be way faster
|
||||
|
||||
stp2.new_process()
|
||||
x2 = stp2()
|
||||
assert np.all(x == x2)
|
||||
|
||||
## STOCPROC FFT ##
|
||||
####################
|
||||
t0 = time.time()
|
||||
stp = sp.StocProc_KLE_tol(r_tau=corr,
|
||||
t_max=t_max,
|
||||
tol=1e-2,
|
||||
ng_fac=4,
|
||||
seed=0)
|
||||
stp = sp.StocProc_FFT(spectral_density, t_max, corr, seed=0)
|
||||
t1 = time.time()
|
||||
dt1 = t1 - t0
|
||||
stp.new_process()
|
||||
x = stp()
|
||||
|
||||
bin_data = pickle.dumps(stp)
|
||||
t0=time.time()
|
||||
stp2 = pickle.loads(bin_data)
|
||||
t1 = time.time()
|
||||
dt2 = t1 - t0
|
||||
assert dt2 / dt1 < 0.1 # loading should be way faster
|
||||
|
||||
stp2.new_process()
|
||||
x2 = stp2()
|
||||
|
||||
assert np.all(x == x2)
|
||||
|
||||
t0 = time.time()
|
||||
stp = sp.StocProc_FFT_tol(spectral_density, t_max, corr, seed=0)
|
||||
t1 = time.time()
|
||||
dt1 = t1-t0
|
||||
|
||||
stp.new_process()
|
||||
x = stp()
|
||||
|
|
@ -222,47 +182,68 @@ def test_stocproc_dump_load():
|
|||
t1 = time.time()
|
||||
dt2 = t1 - t0
|
||||
|
||||
assert dt2/dt1 < 0.1 # loading should be way faster
|
||||
assert dt2 / dt1 < 0.1 # loading should be way faster
|
||||
|
||||
stp2.new_process()
|
||||
x2 = stp2()
|
||||
|
||||
assert np.all(x == x2)
|
||||
|
||||
def test_lorentz_SD(plot=False):
|
||||
_WC_ = 1
|
||||
def lsp(w):
|
||||
return 1/(1 + (w - _WC_)**2)# / np.pi
|
||||
|
||||
def lac(t):
|
||||
return np.exp(- np.abs(t) - 1j*_WC_*t)
|
||||
def test_many(plot=False):
|
||||
t_max = 15
|
||||
num_samples = 5000
|
||||
tol = 5e-2
|
||||
|
||||
sd = spectral_density
|
||||
ac = corr
|
||||
|
||||
stp = sp.StocProc_FFT(sd, t_max, ac, negative_frequencies=False, seed=0, intgr_tol=5e-3, intpl_tol=5e-3)
|
||||
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')
|
||||
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')
|
||||
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')
|
||||
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')
|
||||
stocproc_metatest(stp, num_samples, tol, ac, plot)
|
||||
|
||||
|
||||
t_max = 15
|
||||
num_samples = 15000
|
||||
tol = 3e-2
|
||||
num_samples = 12000
|
||||
tol = 5e-2
|
||||
|
||||
#stp = sp.StocProc_FFT_tol(lsp, t_max, lac, negative_frequencies=True, seed=0, intgr_tol=5e-3, intpl_tol=5e-3)
|
||||
#stocproc_metatest(stp, num_samples, tol, lac, plot)
|
||||
sd = lsd
|
||||
ac = lac
|
||||
|
||||
stp = sp.StocProc_KLE_tol(tol=1e-2,
|
||||
r_tau = lac,
|
||||
t_max = t_max,
|
||||
ng_fac = 2,
|
||||
seed = 0)
|
||||
stocproc_metatest(stp, num_samples, tol, lac, plot)
|
||||
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)
|
||||
|
||||
stp = sp.StocProc_KLE(tol=5e-3, r_tau=ac, t_max=t_max, ng_fac=1, seed=0, diff_method='full', meth='simp')
|
||||
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')
|
||||
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')
|
||||
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')
|
||||
stocproc_metatest(stp, num_samples, tol, ac, plot)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
import logging
|
||||
logging.basicConfig(level=logging.INFO)
|
||||
#logging.basicConfig(level=logging.INFO)
|
||||
# test_stochastic_process_KLE_correlation_function(plot=True)
|
||||
# test_stochastic_process_FFT_correlation_function(plot=True)
|
||||
# test_stochastic_process_KLE_tol_correlation_function(plot=True)
|
||||
# test_stochastic_process_KLE_correlation_function(plot=False)
|
||||
# test_stochastic_process_FFT_correlation_function(plot=False)
|
||||
# test_stocproc_dump_load()
|
||||
|
||||
|
||||
test_lorentz_SD(plot=True)
|
||||
test_many(plot=False)
|
||||
# test_subohmic_SD()
|
||||
pass
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue