diff --git a/doc/source/StocProc_FFT.rst b/doc/source/StocProc_FFT.rst new file mode 100644 index 0000000..4b5fa8d --- /dev/null +++ b/doc/source/StocProc_FFT.rst @@ -0,0 +1,4 @@ +StocProc_FFT_tol +---------------- + +.. autoclass:: stocproc.stocproc.StocProc_FFT_tol \ No newline at end of file diff --git a/doc/source/StocProc_KLE.rst b/doc/source/StocProc_KLE.rst new file mode 100644 index 0000000..d9b55eb --- /dev/null +++ b/doc/source/StocProc_KLE.rst @@ -0,0 +1,9 @@ +StocProc_KLE_tol +---------------- + +.. py:currentmodule:: stocproc + +.. autoclass:: StocProc_KLE_tol + :members: + :inherited-members: + :undoc-members: diff --git a/doc/source/conf.py b/doc/source/conf.py index dbf3001..3a927c1 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -40,6 +40,8 @@ extensions = [ 'sphinx.ext.githubpages', ] +autoclass_content = 'both' + # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] diff --git a/doc/source/example.rst b/doc/source/example.rst new file mode 100644 index 0000000..59fcf1b --- /dev/null +++ b/doc/source/example.rst @@ -0,0 +1,11 @@ +Example +======= + +This example sets up a generator for stochastic processes with exponential +auto correlation. A single process is shown as well as the reconstruction +of the auto correlation using 5000 samples. + +This example is used to generate the figures for the example +section of the :doc:`StocProc Module `. + +.. literalinclude:: ../../examples/example.py \ No newline at end of file diff --git a/doc/source/index.rst b/doc/source/index.rst index 8e25274..ff6f87a 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -1 +1,2 @@ -.. automodule:: stocproc \ No newline at end of file +.. automodule:: stocproc + diff --git a/doc/source/stocproc.rst b/doc/source/stocproc.rst new file mode 100644 index 0000000..33d9ffa --- /dev/null +++ b/doc/source/stocproc.rst @@ -0,0 +1 @@ +.. automodule:: stocproc.stocproc \ No newline at end of file diff --git a/examples/example.py b/examples/example.py new file mode 100644 index 0000000..c20662c --- /dev/null +++ b/examples/example.py @@ -0,0 +1,76 @@ +import sys +import pathlib +# add stocproc package location to path +sys.path.insert(0, pathlib.Path(__file__).parent.parent) +import matplotlib.pyplot as plt +import numpy as np +import stocproc as sp + +_WC_ = 5 +def lsd(w): + """Lorenzian spectral density""" + return 1/(1 + (w - _WC_)**2) + +def lac(t): + """the corresponding Lorenzian correlation function + + note there is a factor of one over pi in the deficition + of the correlation function: + lac(t) = 1/pi int_{-infty}^infty d w lsd(w) exp(-i w t) + """ + return np.exp(- np.abs(t) - 1j*_WC_*t) + +t_max = 10 +print("setup process generator") +stp = sp.StocProc_FFT_tol(lsd, t_max, lac, + negative_frequencies=True, seed=0, + intgr_tol=1e-2, intpl_tol=1e-2) + +fig = plt.figure() +print("generate single process") +stp.new_process() +zt = stp() # get discrete process +t = stp.t # and the natural time axis for the discrete process +plt.plot(stp.t, np.real(stp()), color='k', + label="$\mathrm{real}(z(t))$") +plt.plot(stp.t, np.imag(stp()), color='k', ls='--', + label="$\mathrm{imag}(z(t))$") +plt.xlim([0,10]) +plt.legend(ncol=2, loc='upper right') +plt.title("stochastic process with exponential autocorrelaition") +plt.xlabel("time") +plt.ylabel("process $z(t)$") +plt.grid() +plt.savefig("proc.png") + +ns = 5000 +print("generate {} samples".format(ns)) +# choose time axis 4 time finer than the natural discrete axis +tfine = np.linspace(0, t_max, (stp.num_grid_points - 1) * 4 + 1) +corr = np.zeros(shape=len(tfine), dtype=np.complex128) +# tells that we want to calculate +tref = 2 +# calculates the auto correlation +for i in range(ns): + stp.new_process() + zt = stp(tfine) + corr += zt * np.conj(stp(tref)) +corr /= ns + +fig = plt.figure() +aclab = r"\langle z(t)z^\ast(t_\mathrm{ref})\rangle" +plt.plot(tfine, np.real(lac(tfine-tref)), color='r', + label=r"$\mathrm{{real}}\left({}\right)$".format(aclab)) +plt.plot(tfine, np.imag(lac(tfine-tref)), color='r', ls='--', + label=r"$\mathrm{{imag}}\left({}\right)$".format(aclab)) +plt.plot(tfine, np.real(corr), color='k', + label="exact auto correlation") +plt.plot(tfine, np.imag(corr), color='k', ls='--') +plt.legend(loc='upper right') +plt.title("auto correlation using {} samples".format(ns)) +plt.xlabel("time") +plt.ylabel(r"auto correlation ($t_\mathrm{{ref}}={}$)".format(tref)) +plt.grid() +plt.savefig("ac.png") + +plt.show() diff --git a/stocproc/__init__.py b/stocproc/__init__.py index 3376620..7b055af 100644 --- a/stocproc/__init__.py +++ b/stocproc/__init__.py @@ -1,32 +1,61 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - """ Stochastic Process Module ========================= This module contains two different implementation for generating stochastic processes for a -given auto correlation function. Both methods are based on a time discrete process, however cubic +given auto correlation function (:doc:`Karhunen-Loève expansion ` +and :doc:`Fast-Fourier method `). +Both methods are based on a time discrete process, however cubic spline interpolation is assured to be valid within a given tolerance. -* simulate stochastic processes using Karhunen-Loève expansion :py:func:`stocproc.StocProc_KLE_tol` +.. toctree:: + :maxdepth: 2 - 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. + stocproc + example -* simulate stochastic processes using Fast-Fourier method method :py:func:`stocproc.StocProc_FFT_tol` +Example +------- - 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. +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 lac(t): + # the corresponding Lorenzian correlation function + # note there is a factor of one over pi in the + # deficition of the correlation function: + # lac(t) = 1/pi int_{-infty}^infty d w lsd(w) exp(-i w t) + return np.exp(- np.abs(t) - 1j*_WC_*t) + + t_max = 10 + print("setup process generator") + stp = sp.StocProc_FFT_tol(lsd, t_max, lac, + negative_frequencies=True, seed=0, + 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 `. + +.. 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.0' +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_KLE from .stocproc import StocProc_KLE_tol \ No newline at end of file diff --git a/stocproc/method_fft.py b/stocproc/method_fft.py index e0126b0..e82bbaa 100644 --- a/stocproc/method_fft.py +++ b/stocproc/method_fft.py @@ -1,3 +1,9 @@ +""" + The method_fft module provides convenient function to + setup a stochastic process generator using fft method + + +""" from __future__ import division, print_function from .tools import ComplexInterpolatedUnivariateSpline @@ -10,8 +16,7 @@ from scipy.optimize import brentq from scipy.optimize import basinhopping import sys import warnings - -warnings.simplefilter('error') +#warnings.simplefilter('error') MAX_FLOAT = sys.float_info.max log = logging.getLogger(__name__) @@ -221,7 +226,7 @@ def opt_integral_boundaries(integrand, a, b, t_max, ft_ref, opt_b_only, N): stepsize=0.1*(b-a)) a, b = r.x[0], r.x[1] rd = 10 ** r.fun - log.info("optimization yields max rd {:.3e} and new boundaries [{:.2e},{:.2e}]".format(rd, a, b)) + log.info("optimization with N {} yields max rd {:.3e} and new boundaries [{:.2e},{:.2e}]".format(N, rd, a, b)) return rd, a, b @@ -279,17 +284,19 @@ def calc_ab_N_dx_dt(integrand, intgr_tol, intpl_tol, t_max, a, b, ft_ref, opt_b_ ft_ref = ft_ref, opt_b_only=opt_b_only, N_max = N_max) - dt = get_dt_for_accurate_interpolation(t_max = t_max, - tol = intpl_tol, - ft_ref = ft_ref) + dt_tol = get_dt_for_accurate_interpolation(t_max = t_max, + tol = intpl_tol, + ft_ref = ft_ref) dx = (b-a)/N - N_min = 2*np.pi/dx/dt - if N_min <= N: + dt = 2*np.pi/dx/N + if dt <= dt_tol: log.info("dt criterion fulfilled") return a, b, N, dx, dt else: log.info("down scale dx and dt to match new power of 2 N") + + N_min = 2*np.pi/dx/dt_tol N = 2**int(np.ceil(np.log2(N_min))) diff --git a/stocproc/method_kle.py b/stocproc/method_kle.py index 04f6351..1243e16 100644 --- a/stocproc/method_kle.py +++ b/stocproc/method_kle.py @@ -52,7 +52,7 @@ def solve_hom_fredholm(r, w, eig_val_min): return eig_val, eig_vec -def get_mid_point_weights(t_max, num_grid_points): +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 diff --git a/stocproc/stocproc.py b/stocproc/stocproc.py index d752d58..6d0498a 100644 --- a/stocproc/stocproc.py +++ b/stocproc/stocproc.py @@ -1,6 +1,29 @@ -# -*- coding: utf8 -*- -from __future__ import print_function, division +""" +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_tol `. + +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 +------------------- +simulate stochastic processes using Fast-Fourier method method :py:func:`stocproc.StocProc_FFT_tol` + +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. +""" import numpy as np import time @@ -113,6 +136,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 @@ -134,7 +159,11 @@ class StocProc_KLE(_absStocProc): # this lengthy part will be skipped when init class from dump, as _A and alpha_k will be stored t0 = time.time() - t, w = method_kle.get_simpson_weights_times(t_max, ng_fredholm) + 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: @@ -229,11 +258,49 @@ class StocProc_KLE(_absStocProc): 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 + + .. math:: Z(t) = \sum_i \sqrt{\lambda_i} Y_i u_i(t) + + where :math:`Y_i` and independent complex valued Gaussian random variables with variance one + (:math:`\langle Y_i Y_j \rangle = \delta_{ij}`) and :math:`\lambda_i`, :math:`u_i(t)` are + eigenvalues / eigenfunctions of the following Fredholm equation + + .. math:: \int_0^{t_\mathrm{max}} \mathrm{d}s R(t-s) u_i(s) = \lambda_i u_i(t) + + 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 + + + - 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 + same as StocProc_KLE except that ng_fredholm is determined from given tolerance + + bla bla + """ - def __init__(self, tol=1e-2, **kwargs): + def __init__(self, r_tau, t_max, tol=1e-2, ng_fac=4, seed=None, k=3, align_eig_vec=False): + """this is init + + :param r_tau: + :param t_max: + :param tol: + :param ng_fac: + :param seed: + :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) @@ -339,6 +406,7 @@ class StocProc_FFT_tol(_absStocProc): N_max = 2**24) log.info("required tol result in N {}".format(N)) + assert abs(2*np.pi - N*dx*dt) < 1e-12 num_grid_points = int(np.ceil(t_max/dt))+1 t_max = (num_grid_points-1)*dt diff --git a/stocproc/tools.py b/stocproc/tools.py index 5f918f3..e738c9b 100644 --- a/stocproc/tools.py +++ b/stocproc/tools.py @@ -20,13 +20,22 @@ stocproc_key_type = namedtuple(typename = 'stocproc_key_type', class ComplexInterpolatedUnivariateSpline(object): - def __init__(self, x, y, k=2): + def __init__(self, x, y, k=3): self.re_spline = InterpolatedUnivariateSpline(x, np.real(y)) self.im_spline = InterpolatedUnivariateSpline(x, np.imag(y)) def __call__(self, t): return self.re_spline(t) + 1j * self.im_spline(t) + +def complex_quad(func, a, b, **kw_args): + func_re = lambda t: np.real(func(t)) + func_im = lambda t: np.imag(func(t)) + I_re = quad(func_re, a, b, **kw_args)[0] + I_im = quad(func_im, a, b, **kw_args)[0] + + return I_re + 1j * I_im + def auto_correlation_numpy(x, verbose=1): warn("use 'auto_correlation' instead", DeprecationWarning) diff --git a/tests/test_method_kle.py b/tests/test_method_kle.py new file mode 100644 index 0000000..335b3d5 --- /dev/null +++ b/tests/test_method_kle.py @@ -0,0 +1,196 @@ +import sys +import os + +import numpy as np +import math +from scipy.special import gamma as gamma_func +import scipy.integrate as sp_int +try: + import matplotlib.pyplot as plt +except ImportError: + print("matplotlib not found -> any plotting will crash") + + +import pathlib +p = pathlib.PosixPath(os.path.abspath(__file__)) +sys.path.insert(0, str(p.parent.parent)) + +import stocproc as sp +from stocproc import tools +from stocproc import stocproc_c +from scipy.integrate import quad +import logging + +def test_solve_fredholm(): + _WC_ = 2 + def lac(t): + return np.exp(- np.abs(t) - 1j*_WC_*t) + + t_max = 10 + + for ng in range(11, 450, 30): + t, w = sp.method_kle.get_mid_point_weights_times(t_max, ng) + r = lac(t.reshape(-1,1)-t.reshape(1,-1)) + _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0) + + ng_fac = 4 + ng_fine = ng_fac * (ng - 1) + 1 + tfine = np.linspace(0, t_max, ng_fine) + bcf_n_plus = lac(tfine - tfine[0]) + alpha_k = np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus)) + + u_i_all_t = stocproc_c.eig_func_all_interp(delta_t_fac=ng_fac, + time_axis=t, + alpha_k=alpha_k, + weights=w, + eigen_val=_eig_val, + eigen_vec=_eig_vec) + + u_i_all_ast_s = np.conj(u_i_all_t) # (N_gp, N_ev) + num_ev = len(_eig_val) + tmp = _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] = alpha_k[idx:idx + ng_fine] + + err = np.max(np.abs(recs_bcf - refc_bcf) / np.abs(refc_bcf)) + plt.plot(ng, err, marker='o', color='r') + + ng += 2 + t, w = sp.method_kle.get_simpson_weights_times(t_max, ng) + r = lac(t.reshape(-1, 1) - t.reshape(1, -1)) + _eig_val2, _eig_vec2 = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0) + #print(np.max(np.abs(_eig_val - _eig_val2)), np.max(np.abs(_eig_vec - _eig_vec2))) + _eig_val, _eig_vec = _eig_val2, _eig_vec2 + + ng_fac = 4 + ng_fine = ng_fac * (ng - 1) + 1 + tfine = np.linspace(0, t_max, ng_fine) + bcf_n_plus = lac(tfine - tfine[0]) + alpha_k = np.hstack((np.conj(bcf_n_plus[-1:0:-1]), bcf_n_plus)) + + u_i_all_t2 = stocproc_c.eig_func_all_interp(delta_t_fac=ng_fac, + time_axis=t, + alpha_k=alpha_k, + weights=w, + eigen_val=_eig_val, + eigen_vec=_eig_vec) + #print(np.max(np.abs(u_i_all_t - u_i_all_t2))) + u_i_all_t = u_i_all_t2 + #print() + + + u_i_all_ast_s = np.conj(u_i_all_t) # (N_gp, N_ev) + num_ev = len(_eig_val) + tmp = _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] = alpha_k[idx:idx + ng_fine] + + err2 = np.max(np.abs(recs_bcf - refc_bcf) / np.abs(refc_bcf)) + plt.plot(ng, err2, marker='o', color='k') + print(err, err2) + + plt.yscale('log') + plt.grid() + plt.show() + +def test_solve_fredholm_reconstr_ac(): + """ + here we see that the reconstruction quality is independent of the integration weights + + differences occur when checking validity of the interpolated time continuous Fredholm equation + """ + _WC_ = 2 + def lac(t): + return np.exp(- np.abs(t) - 1j*_WC_*t) + t_max = 10 + for ng in range(11,500,30): + print(ng) + t, w = sp.method_kle.get_mid_point_weights_times(t_max, ng) + r = lac(t.reshape(-1,1)-t.reshape(1,-1)) + _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0) + _eig_vec_ast = np.conj(_eig_vec) # (N_gp, N_ev) + tmp = _eig_val.reshape(1, -1) * _eig_vec # (N_gp, N_ev) + recs_bcf = np.tensordot(tmp, _eig_vec_ast, axes=([1], [1])) + rd = np.max(np.abs(recs_bcf - r) / np.abs(r)) + assert rd < 1e-10 + + + t, w = sp.method_kle.get_simpson_weights_times(t_max, ng) + r = lac(t.reshape(-1, 1) - t.reshape(1, -1)) + _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0) + _eig_vec_ast = np.conj(_eig_vec) # (N_gp, N_ev) + tmp = _eig_val.reshape(1, -1) * _eig_vec # (N_gp, N_ev) + recs_bcf = np.tensordot(tmp, _eig_vec_ast, axes=([1], [1])) + rd = np.max(np.abs(recs_bcf - r) / np.abs(r)) + assert rd < 1e-10 + +def test_solve_fredholm_interp_eigenfunc(plot=False): + """ + here we take the discrete eigenfunctions of the Fredholm problem + and use qubic interpolation to check the integral equality. + + the difference between the midpoint weights and simpson weights become + visible. Although the simpson integration yields on average a better performance + there are high fluctuation in the error. + """ + _WC_ = 2 + def lac(t): + return np.exp(- np.abs(t) - 1j*_WC_*t) + + t_max = 10 + ng = 81 + ngfac = 2 + tfine = np.linspace(0, t_max, (ng-1)*ngfac+1) + + if plot: + fig, ax = plt.subplots(nrows=2, ncols=2, sharey=True, sharex=True) + ax = ax.flatten() + + for idx in range(4): + t, w = sp.method_kle.get_mid_point_weights_times(t_max, ng) + r = lac(t.reshape(-1, 1) - t.reshape(1, -1)) + _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0) + u0 = tools.ComplexInterpolatedUnivariateSpline(t, _eig_vec[:,idx]) + lam0 = _eig_val[idx] + err_mp = [] + for ti in tfine: + I = tools.complex_quad(lambda s: lac(ti-s)*u0(s), 0, t_max, limit=500) + err_mp.append(np.abs(I - lam0*u0(ti))) + if plot: + axc = ax[idx] + axc.plot(tfine, err_mp, color='r', label='midp') + + t, w = sp.method_kle.get_simpson_weights_times(t_max, ng) + r = lac(t.reshape(-1, 1) - t.reshape(1, -1)) + _eig_val, _eig_vec = sp.method_kle.solve_hom_fredholm(r, w, eig_val_min=0) + u0 = tools.ComplexInterpolatedUnivariateSpline(t, _eig_vec[:, idx]) + lam0 = _eig_val[idx] + err_sp = [] + for ti in tfine: + I = tools.complex_quad(lambda s: lac(ti - s) * u0(s), 0, t_max, limit=500) + err_sp.append(np.abs(I - lam0 * u0(ti))) + if plot: + axc.plot(tfine, err_sp, color='k', label='simp') + axc.set_yscale('log') + axc.plot(tfine[::ngfac], err_sp[::ngfac], ls='', marker='x', color='k') + axc.set_title("eigen function # {}".format(idx)) + + assert max(err_mp) > max(err_sp) + + if plot: + fig.suptitle("np.abs(int R(t-s)u_i(s) - lam_i * u_i(t))") + plt.show() + + +if __name__ == "__main__": + # test_solve_fredholm() + # test_solve_fredholm_reconstr_ac() + test_solve_fredholm_interp_eigenfunc(plot=True) \ No newline at end of file diff --git a/tests/test_stocproc.py b/tests/test_stocproc.py index 9b8f21d..5e6c402 100644 --- a/tests/test_stocproc.py +++ b/tests/test_stocproc.py @@ -40,15 +40,16 @@ def spectral_density(omega): def stocproc_metatest(stp, num_samples, tol, corr, plot): print("generate samples") - t = np.linspace(0, stp.t_max, 487) - x_t_array_KLE = np.empty(shape=(num_samples, 487), dtype=np.complex128) + N = 287 + t = np.linspace(0, stp.t_max, N) + x_t_array_KLE = np.empty(shape=(num_samples, N), dtype=np.complex128) for i in range(num_samples): stp.new_process() x_t_array_KLE[i] = stp(t) autoCorr_KLE_conj, autoCorr_KLE_not_conj = sp.tools.auto_correlation(x_t_array_KLE) - ac_true = corr(t.reshape(487, 1) - t.reshape(1, 487)) + ac_true = corr(t.reshape(N, 1) - t.reshape(1, N)) max_diff_conj = np.max(np.abs(ac_true - autoCorr_KLE_conj)) print("max diff : {:.2e}".format(max_diff_conj)) @@ -63,32 +64,32 @@ def stocproc_metatest(stp, num_samples, tol, corr, plot): v_min_imag = np.floor(np.min(np.imag(ac_true))) v_max_imag = np.ceil(np.max(np.imag(ac_true))) - fig, ax = plt.subplots(nrows=4, ncols=2, figsize=(10, 14)) + fig, ax = plt.subplots(nrows=2, ncols=4, figsize=(14, 10)) ax[0, 0].set_title(r"exact $\mathrm{re}\left(\langle x(t) x^\ast(s) \rangle\right)$") ax[0, 0].imshow(np.real(ac_true), interpolation='none', vmin=v_min_real, vmax=v_max_real, cmap="seismic") - ax[0, 1].set_title(r"exact $\mathrm{im}\left(\langle x(t) x^\ast(s) \rangle\right)$") - ax[0, 1].imshow(np.imag(ac_true), interpolation='none', vmin=v_min_imag, vmax=v_max_imag, cmap="seismic") + ax[1, 0].set_title(r"exact $\mathrm{im}\left(\langle x(t) x^\ast(s) \rangle\right)$") + ax[1, 0].imshow(np.imag(ac_true), interpolation='none', vmin=v_min_imag, vmax=v_max_imag, cmap="seismic") - ax[1, 0].set_title(r"KLE $\mathrm{re}\left(\langle x(t) x^\ast(s) \rangle\right)$") - ax[1, 0].imshow(np.real(autoCorr_KLE_conj), interpolation='none', vmin=v_min_real, vmax=v_max_real, + ax[0, 1].set_title(r"$\mathrm{re}\left(\langle x(t) x^\ast(s) \rangle\right)$") + ax[0, 1].imshow(np.real(autoCorr_KLE_conj), interpolation='none', vmin=v_min_real, vmax=v_max_real, cmap="seismic") - ax[1, 1].set_title(r"KLE $\mathrm{im}\left(\langle x(t) x^\ast(s) \rangle\right)$") + ax[1, 1].set_title(r"$\mathrm{im}\left(\langle x(t) x^\ast(s) \rangle\right)$") ax[1, 1].imshow(np.imag(autoCorr_KLE_conj), interpolation='none', vmin=v_min_imag, vmax=v_max_imag, cmap="seismic") - ax[2, 0].set_title(r"KLE $\mathrm{re}\left(\langle x(t) x(s) \rangle\right)$") - ax[2, 0].imshow(np.real(autoCorr_KLE_not_conj), interpolation='none', vmin=v_min_real, vmax=v_max_real, + ax[0, 2].set_title(r"$\mathrm{re}\left(\langle x(t) x(s) \rangle\right)$") + ax[0, 2].imshow(np.real(autoCorr_KLE_not_conj), interpolation='none', vmin=v_min_real, vmax=v_max_real, cmap="seismic") - ax[2, 1].set_title(r"KLE $\mathrm{im}\left(\langle x(t) x(s) \rangle\right)$") - ax[2, 1].imshow(np.imag(autoCorr_KLE_not_conj), interpolation='none', vmin=v_min_imag, vmax=v_max_imag, + ax[1, 2].set_title(r"$\mathrm{im}\left(\langle x(t) x(s) \rangle\right)$") + ax[1, 2].imshow(np.imag(autoCorr_KLE_not_conj), interpolation='none', vmin=v_min_imag, vmax=v_max_imag, cmap="seismic") - ax[3, 0].set_title(r"abs diff $\langle x(t) x^\ast(s) \rangle$") - cax = ax[3, 0].imshow(np.log10(np.abs(autoCorr_KLE_conj - ac_true)), interpolation='none', cmap="inferno") - fig.colorbar(cax, ax=ax[3, 0]) - ax[3, 1].set_title(r"abs diff $\langle x(t) x(s) \rangle$") - cax = ax[3, 1].imshow(np.log10(np.abs(autoCorr_KLE_not_conj)), interpolation='none', cmap="inferno") - fig.colorbar(cax, ax=ax[3, 1]) + ax[0, 3].set_title(r"abs diff $\langle x(t) x^\ast(s) \rangle$") + cax = ax[0, 3].imshow(np.log10(np.abs(autoCorr_KLE_conj - ac_true)), interpolation='none', cmap="inferno") + fig.colorbar(cax, ax=ax[0, 3]) + ax[1, 3].set_title(r"abs diff $\langle x(t) x(s) \rangle$") + cax = ax[1, 3].imshow(np.log10(np.abs(autoCorr_KLE_not_conj)), interpolation='none', cmap="inferno") + fig.colorbar(cax, ax=ax[1, 3]) plt.tight_layout() plt.show() @@ -238,126 +239,18 @@ def test_lorentz_SD(plot=False): t_max = 15 - num_samples = 5000 + num_samples = 15000 tol = 3e-2 - stp = sp.StocProc_KLE(r_tau=lac, - t_max=t_max, - ng_fredholm=1025, - ng_fac=2, - seed=0) - t = stp.t - stp.new_process() - plt.plot(t, np.abs(stp())) - - - stp = sp.StocProc_FFT_tol(lsp, t_max, lac, negative_frequencies=True, seed=0, intgr_tol=1e-2, intpl_tol=1e-9) - t = stp.t - stp.new_process() - plt.plot(t, np.abs(stp()), ls='', marker='.') - - stp = sp.StocProc_FFT_tol(lsp, t_max, lac, negative_frequencies=True, seed=0, intgr_tol=1e-2, intpl_tol=1e-11) - t = stp.t - stp.new_process() - plt.plot(t, np.abs(stp()), ls='', marker='.') - - - plt.show() - return - - print("generate samples") - t = np.linspace(0, stp.t_max, 487) - t_fine = np.linspace(0, stp.t_max, 4870) - - t_ref_list = [0,1,7] - - - - # for num_samples in [100, 1000]: - # x_t0 = np.zeros(shape=(3), dtype=np.complex128) - # for i in range(num_samples): - # stp.new_process() - # x_t = stp(t) - # for j, ti in enumerate(t_ref_list): - # x_t0[j] += np.conj(stp(ti)) * stp(ti) - # x_t0 /= num_samples - # print(np.abs(1-x_t0)) - - num_samples = 1000 - x_t_array = np.zeros(shape=(487, 3), dtype=np.complex128) - for i in range(num_samples): - stp.new_process() - x_t = stp(t) - for j, ti in enumerate(t_ref_list): - x_t_array[:, j] += x_t * np.conj(stp(ti)) - - x_t_array /= num_samples - for j, ti in enumerate(t_ref_list): - p, = plt.plot(t, np.abs(x_t_array[:, j])) - plt.plot(t_fine, np.abs(lac(t_fine-ti)), color=p.get_color()) - - plt.show() - - - - - #num_samples = 1000 - #stp = sp.StocProc_FFT_tol(lsp, t_max, lac, negative_frequencies=True, seed=0, intgr_tol=1e-2, intpl_tol=1e-2) - #stocproc_metatest(stp, num_samples, tol, lambda tau:lac(tau)/np.pi, plot) - - - -def test_subohmic_SD(plot=False): - - def lsp(omega): - return omega ** _S_ * np.exp(-omega) - - def lac(tau): - return (1 + 1j*(tau))**(-(_S_+1)) * _GAMMA_S_PLUS_1 / np.pi - - - t_max = 15 - stp = sp.StocProc_KLE(r_tau=lac, - t_max=t_max, - ng_fredholm=1025, - ng_fac=2, - seed=0) - t = stp.t - stp.new_process() - plt.plot(t, np.abs(stp())) - - - stp = sp.StocProc_FFT_tol(lsp, t_max, lac, negative_frequencies=False, seed=0, intgr_tol=1e-3, intpl_tol=1e-3) - t = stp.t - stp.new_process() - plt.plot(t, np.abs(stp())) - plt.show() - return - - - - print("generate samples") - t = np.linspace(0, stp.t_max, 487) - t_fine = np.linspace(0, stp.t_max, 4870) - - t_ref_list = [0,1,7] - - num_samples = 1000 - x_t_array = np.zeros(shape=(487, 3), dtype=np.complex128) - for i in range(num_samples): - stp.new_process() - x_t = stp(t) - for j, ti in enumerate(t_ref_list): - x_t_array[:, j] += x_t * np.conj(stp(ti)) - - x_t_array /= num_samples - for j, ti in enumerate(t_ref_list): - p, = plt.plot(t, np.abs(x_t_array[:, j])) - plt.plot(t_fine, np.abs(lac(t_fine-ti)), color=p.get_color()) - - plt.show() - + #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) + 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) if __name__ == "__main__":