From 8fdf4ef77785f05d6914de179464c2feb08c253d Mon Sep 17 00:00:00 2001 From: cimatosa Date: Thu, 18 Jun 2015 10:30:22 +0200 Subject: [PATCH] intermed work --- __init__.py | 1 + stocproc.py | 256 ++++++++++++++++++++++++++++++++--------------- test_gquad.py | 7 +- test_stocproc.py | 137 +++++++++++++++++++------ 4 files changed, 292 insertions(+), 109 deletions(-) diff --git a/__init__.py b/__init__.py index e69de29..1fc35d4 100644 --- a/__init__.py +++ b/__init__.py @@ -0,0 +1 @@ +from .stocproc import * \ No newline at end of file diff --git a/stocproc.py b/stocproc.py index 08bc2a8..bc1e2c7 100644 --- a/stocproc.py +++ b/stocproc.py @@ -182,18 +182,15 @@ class StocProc(object): # solve discrete Fredholm equation # eig_val = lambda # eig_vec = u(t) - self._eig_val, self._eig_vec = solve_hom_fredholm(r, w, sig_min**2) - self._sqrt_eig_val = np.sqrt(self._eig_val) - self._num_ev = len(self._eig_val) - self._A = self._w.reshape(self._num_gp,1) * self._eig_vec / self._sqrt_eig_val.reshape(1, self._num_ev) - + self._eig_val, self._eig_vec = solve_hom_fredholm(r, w, sig_min**2, verbose=self.verbose) else: self.__load(fname) - self.__calc_missing() + + self.__calc_missing() self.my_cache_decorator = functools.lru_cache(maxsize=cache_size, typed=False) self.x = self.my_cache_decorator(self._x) - self.new_process(seed) + self.new_process(seed = seed) @classmethod def new_instance_by_name(cls, name, r_tau, t_max, ng, seed, sig_min): @@ -211,19 +208,19 @@ class StocProc(object): ob.name = name return ob @classmethod - def new_instance_with_trapezoidal_weights(cls, r_tau, t_max, ng, seed, sig_min): + def new_instance_with_trapezoidal_weights(cls, r_tau, t_max, ng, seed, sig_min, verbose=1): t, w = get_trapezoidal_weights_times(t_max, ng) - return cls(r_tau, t, w, seed, sig_min) + return cls(r_tau, t, w, seed, sig_min, verbose=verbose) @classmethod - def new_instance_with_mid_point_weights(cls, r_tau, t_max, ng, seed, sig_min): + def new_instance_with_mid_point_weights(cls, r_tau, t_max, ng, seed, sig_min, verbose=1): t, w = get_mid_point_weights(t_max, ng) - return cls(r_tau, t, w, seed, sig_min) + return cls(r_tau, t, w, seed, sig_min, verbose=verbose) @classmethod - def new_instance_with_gauss_legendre_weights(cls, r_tau, t_max, ng, seed, sig_min): + def new_instance_with_gauss_legendre_weights(cls, r_tau, t_max, ng, seed, sig_min, verbose=1): t, w = gquad.gauss_nodes_weights_legendre(n=ng, low=0, high=t_max) - return cls(r_tau, t, w, seed, sig_min) + return cls(r_tau, t, w, seed, sig_min, verbose=verbose) def __load(self, fname): with open(fname, 'rb') as f: @@ -234,7 +231,8 @@ class StocProc(object): self._num_gp = len(self._s) self._sqrt_eig_val = np.sqrt(self._eig_val) self._num_ev = len(self._eig_val) - self._A = self._w.reshape(self._num_gp,1) * self._eig_vec / self._sqrt_eig_val.reshape(1, self._num_ev) + self._w /= np.sqrt(2) + self._A = self._w.reshape(self._num_gp,1) * self._eig_vec / self._sqrt_eig_val.reshape(1, self._num_ev) def __dump(self, fname): @@ -259,7 +257,7 @@ class StocProc(object): else: return 'unknown' - def new_process(self, seed = None): + def new_process(self, yi=None, seed=None): r"""setup new process Generates a new set of independent normal random variables :math:`Y_i` @@ -274,7 +272,14 @@ class StocProc(object): np.random.seed(seed) self.clear_cache() - self._Y = 1/np.sqrt(2)*np.random.normal(size=2*self._num_ev).view(np.complex).reshape(self._num_ev,1) + if yi is None: + if self.verbose > 1: + print("generate samples ...") + self._Y = np.random.normal(size=2*self._num_ev).view(np.complex).reshape(self._num_ev,1) + if self.verbose > 1: + print("done!") + else: + self._Y = yi.reshape(self._num_ev,1) def x_for_initial_time_grid(self): @@ -285,8 +290,14 @@ class StocProc(object): equation. This is equivalent to calling :py:func:`stochastic_process_kle` with the same weights :math:`w_i` and time grid points :math:`s_i`. """ - tmp = self._Y * self._sqrt_eig_val.reshape(self._num_ev,1) - return np.tensordot(tmp, self._eig_vec, axes=([0],[1])).flatten() + tmp = self._Y / np.sqrt(2) * self._sqrt_eig_val.reshape(self._num_ev,1) + if self.verbose > 1: + print("calc process via matrix prod ...") + res = np.tensordot(tmp, self._eig_vec, axes=([0],[1])).flatten() + if self.verbose > 1: + print("done!") + + return res def time_grid(self): return self._s @@ -295,9 +306,8 @@ class StocProc(object): if isinstance(t, np.ndarray): return self.x_t_array(t) else: - return self._x(t) + return self.x(t) - @functools.lru_cache(maxsize=1024, typed=False) def _x(self, t): # _Y # (N_ev, 1 ) tmp = self._Y*self._r_tau(t-self._s.reshape(1, self._num_gp)) @@ -319,7 +329,12 @@ class StocProc(object): # (N_ev, N_gp, N_t) # A # (N_gp, N_ev) # A_j,i = w_j / sqrt(lambda_i) u_i(s_j) - return np.tensordot(tmp, self._A, axes=([1,0],[0,1])) + if self.verbose > 1: + print("calc process via matrix prod ...") + res = np.tensordot(tmp, self._A, axes=([1,0],[0,1])) + if self.verbose > 1: + print("done!") + return res def u_i(self, t_array, i): @@ -340,7 +355,7 @@ class StocProc(object): # (N_gp, N_t) # A # (N_gp, N_ev) # A_j,i = w_j / sqrt(lambda_i) u_i(s_j) - return 1/self._sqrt_eig_val[i]*np.tensordot(tmp, self._A[:,i], axes=([0],[0])) + return np.sqrt(2)/self._sqrt_eig_val[i]*np.tensordot(tmp, self._A[:,i], axes=([0],[0])) def u_i_all(self, t_array): r"""get all eigenfunctions @@ -362,7 +377,7 @@ class StocProc(object): # (N_gp, N_t) # A # (N_gp, N_ev) # A_j,i = w_j / sqrt(lambda_i) u_i(s_j) - return np.tensordot(tmp, 1/self._sqrt_eig_val.reshape(1,self._num_ev) * self._A, axes=([0],[0])) + return np.tensordot(tmp, np.sqrt(2)/self._sqrt_eig_val.reshape(1,self._num_ev) * self._A, axes=([0],[0])) def eigen_vector_i(self, i): r"""Returns the i-th eigenvector (solution of the discrete Fredhom equation)""" @@ -425,19 +440,83 @@ class StocProc(object): tmp = lambda_i_all.reshape(1, self._num_ev) * u_i_all_t #(N_gp, N_ev) return np.tensordot(tmp, u_i_all_ast_s, axes=([1],[1]))[:,0] +class StocProc_KLE(StocProc): + def __init__(self, r_tau, t_max, num_grid_points, ng_fredholm=None, seed=None, sig_min=1e-5, verbose=1): + """ + class to simulate stocastic processes using KLE method + - solve fredholm equation on grid with ng_fregholm nodes (trapezoidal_weights) + - calculate process (using interpolation solution of fredholm equation) with num_grid_points nodes + - ivoke cubic spline interpolator when calling + """ + + if ng_fredholm is None: + ng_fredholm = num_grid_points + + if num_grid_points < ng_fredholm: + print("WARNING: found 'num_grid_points < ng_fredholm' -> set 'num_grid_points == ng_fredholm'") + + if ng_fredholm == num_grid_points: + self.kle_interp = False + else: + self.kle_interp = True + + self.stocproc = StocProc.new_instance_with_trapezoidal_weights(r_tau = r_tau, + t_max = t_max, + ng = ng_fredholm, + sig_min = sig_min, + seed = seed, + verbose = verbose) + self._r_tau = self.stocproc._r_tau + self._s = self.stocproc._s + self._A = self.stocproc._A + self.stocproc.verbose = verbose + self.verbose = verbose + self.t = np.linspace(0, t_max, num_grid_points) + self._z = self.__cal_z() + if self.verbose > 1: + print("setup interpolator ...") + self.interpolator = ComplexInterpolatedUnivariateSpline(self.t, self._z, k=3) + if self.verbose > 1: + print("done!") + + def __cal_z(self): + if self.kle_interp: + return self.stocproc.x_t_array(self.t) + else: + return self.stocproc.x_for_initial_time_grid() + + def __call__(self, t): + return self.interpolator(t) + + def get_time(self): + return self.t + + def get_z(self): + return self._z + + def new_process(self, yi = None, seed = None): + self.stocproc.new_process(yi=yi, seed=seed) + self._z = self.__cal_z() + if self.verbose > 1: + print("setup interpolator ...") + self.interpolator = ComplexInterpolatedUnivariateSpline(self.t, self._z, k=3) + if self.verbose > 1: + print("done!") + + # def reconst_corr_zero(self, t_array): class StocProc_FFT(object): r""" Simulate Stochastic Process using FFT method """ - def __init__(self, spectral_density, t_max, num_grid_points, seed = None, verbose=1): - self.verbose = verbose - self.t_max = t_max + def __init__(self, spectral_density, t_max, num_grid_points, seed=None, verbose=1): + self.verbose = verbose + self.t_max = t_max self.num_grid_points = num_grid_points - self.n_dft = self.num_grid_points * 2 - 1 - delta_t = self.t_max / (self.num_grid_points-1) - self.delta_omega = 2 * np.pi / (delta_t * self.n_dft) + self.n_dft = self.num_grid_points * 2 - 1 + delta_t = self.t_max / (self.num_grid_points-1) + self.delta_omega = 2 * np.pi / (delta_t * self.n_dft) #omega axis omega = self.delta_omega*np.arange(self.n_dft) @@ -445,36 +524,56 @@ class StocProc_FFT(object): self.sqrt_spectral_density_times_sqrt_delta_omega_over_sqrt_2 = np.sqrt(spectral_density(omega)) * np.sqrt(self.delta_omega) / np.sqrt(2) if self.verbose > 0: - print(" omega_max : {:.2}".format(self.delta_omega * self.n_dft)) - print(" delta_omega: {:.2}".format(self.delta_omega)) + print("stoc proc fft, spectral density sampling information:") + print(" t_max :", (t_max)) + print(" ng :", (num_grid_points)) + print(" omega_max :", (self.delta_omega * self.n_dft)) + print(" delta_omega:", (self.delta_omega)) + + self.interpolator = None + self.t = np.linspace(0, self.t_max, self.num_grid_points) self.new_process(seed = seed) + + def __call__(self, t): + if self.interpolator is None: + if self.verbose > 1: + print("setup interpolator ...") + self.interpolator = ComplexInterpolatedUnivariateSpline(self.t, self._z, k=3) + if self.verbose > 1: + print("done!") + return self.interpolator(t) + + def get_time(self): + return self.t + + def get_z(self): + return self._z - def new_process(self, seed = None): + def new_process(self, yi = None, seed = None): self.interpolator = None if seed != None: if self.verbose > 0: print("use seed", seed) np.random.seed(seed) - #random complex normal samples - yi = np.random.normal(size = 2*self.n_dft).view(np.complex) + if yi is None: + #random complex normal samples + if self.verbose > 1: + print("generate samples ...") + yi = np.random.normal(size = 2*self.n_dft).view(np.complex) + if self.verbose > 1: + print("done") #each row contain a different integrand weighted_integrand = self.sqrt_spectral_density_times_sqrt_delta_omega_over_sqrt_2 * yi #compute integral using fft routine - self.z_ast = np.fft.fft(weighted_integrand)[0:self.num_grid_points] + if self.verbose > 1: + print("calc process via fft ...") + self._z = np.fft.fft(weighted_integrand)[0:self.num_grid_points] + if self.verbose > 1: + print("done") + - def x_for_initial_time_grid(self): - return self.z_ast - - def __call__(self, t): - if self.interpolator is None: - self.interpolator = ComplexInterpolatedUnivariateSpline(self.get_time_axis(), self.z_ast, k=3) - return self.interpolator(t) - - def get_time_axis(self): - #corresponding time axis - return np.linspace(0, self.t_max, self.num_grid_points) def _mean_error(r_t_s, r_t_s_exact): @@ -484,64 +583,65 @@ def _mean_error(r_t_s, r_t_s_exact): :return: the mean error ``err`` """ - len_t, len_s = r_t_s.shape - abs_sqare = np.abs(r_t_s - r_t_s_exact)**2 - abs_sqare[0,:] /= 2 - abs_sqare[-1,:] /= 2 - - err = np.sum(abs_sqare, axis = 0) / len_t + err = np.mean(np.abs(r_t_s - r_t_s_exact), axis = 0) return err def _max_error(r_t_s, r_t_s_exact): - return np.max(np.abs(r_t_s - r_t_s_exact)) + return np.max(np.abs(r_t_s - r_t_s_exact), axis = 0) def _max_rel_error(r_t_s, r_t_s_exact): return np.max(np.abs(r_t_s - r_t_s_exact) / np.abs(r_t_s_exact)) -def auto_grid_points(r_tau, t_max, ng_interpolation, tol = 1e-8, err_method = _max_error, name = 'mid_point', sig_min = 1e-4): +def auto_grid_points(r_tau, t_max, tol = 1e-8, err_method = _max_error, name = 'mid_point', sig_min = 1e-4): err = 1 - ng = 1 + ng = 64 seed = None - t_large = np.linspace(0, t_max, ng_interpolation) + err_method_name = err_method.__name__ print("start auto_grid_points, determine ng ...") #exponential increase to get below error threshold while err > tol: ng *= 2 + ng_fine = ng*10 + t_fine = np.linspace(0, t_max, ng_fine) + print("#"*40) + print("new process with {} weights ...".format(name)) stoc_proc = StocProc.new_instance_by_name(name, r_tau, t_max, ng, seed, sig_min) - print(" new process with {} weights".format(stoc_proc.get_name())) - r_t_s = stoc_proc.recons_corr(t_large) - r_t_s_exact = r_tau(t_large.reshape(ng_interpolation,1) - t_large.reshape(1, ng_interpolation)) - - err = err_method(r_t_s, r_t_s_exact) - - print(" ng {} -> err {:.2e}".format(ng, err)) + print("reconstruct correlation function ({} points)...".format(ng_fine)) + r_t_s = stoc_proc.recons_corr(t_fine) + print("calculate exact correlation function ...") + r_t_s_exact = r_tau(t_fine.reshape(ng_fine,1) - t_fine.reshape(1, ng_fine)) + print("calculate error using {} ...".format(err_method_name)) + err = np.max(err_method(r_t_s, r_t_s_exact)) + print("ng {} -> err {:.3e}".format(ng, err)) ng_low = ng // 2 ng_high = ng while (ng_high - ng_low) > 1: print("#"*40) - print(" ng_l", ng_low) - print(" ng_h", ng_high) + print("ng_l", ng_low) + print("ng_h", ng_high) ng = (ng_low + ng_high) // 2 - print(" ng", ng) - + print("ng", ng) + ng_fine = ng*10 + t_fine = np.linspace(0, t_max, ng_fine) + print("new process with {} weights ...".format(name)) stoc_proc = StocProc.new_instance_by_name(name, r_tau, t_max, ng, seed, sig_min) - - r_t_s = stoc_proc.recons_corr(t_large) - r_t_s_exact = r_tau(t_large.reshape(ng_interpolation,1) - t_large.reshape(1, ng_interpolation)) - - err = err_method(r_t_s, r_t_s_exact) - - print(" ng {} -> err {:.2e}".format(ng, err)) + print("reconstruct correlation function ({} points)...".format(ng_fine)) + r_t_s = stoc_proc.recons_corr(t_fine) + print("calculate exact correlation function ...") + r_t_s_exact = r_tau(t_fine.reshape(ng_fine,1) - t_fine.reshape(1, ng_fine)) + print("calculate error using {} ...".format(err_method_name)) + err = np.max(err_method(r_t_s, r_t_s_exact)) + print("ng {} -> err {:.3e}".format(ng, err)) if err > tol: - print(" err > tol!") - print(" ng_l -> ", ng) + print(" err > tol!") + print(" ng_l -> ", ng) ng_low = ng else: - print(" err <= tol!") - print(" ng_h -> ", ng) + print(" err <= tol!") + print(" ng_h -> ", ng) ng_high = ng diff --git a/test_gquad.py b/test_gquad.py index ae2a7a1..9739594 100644 --- a/test_gquad.py +++ b/test_gquad.py @@ -1,9 +1,14 @@ -import gquad from scipy.integrate import quad import numpy as np import numpy.polynomial as pln import pytest +import os +import sys +path = os.path.dirname(__file__) +sys.path.append(path) +import gquad + def scp_laguerre(p1,p2): f = lambda x: p1(x)*p2(x)*np.exp(-x) diff --git a/test_stocproc.py b/test_stocproc.py index 9377b9b..24a9a9b 100644 --- a/test_stocproc.py +++ b/test_stocproc.py @@ -182,6 +182,7 @@ def test_func_vs_class_KLE_FFT(): sig_min = sig_min) x_t_array_class = stoc_proc.x_for_initial_time_grid() + print(np.max(np.abs(x_t_array_func - x_t_array_class))) assert np.all(x_t_array_func == x_t_array_class), "stochastic_process_kle vs. StocProc Class not identical" x_t_array_func, t = sp.stochastic_process_fft(spectral_density = J, @@ -195,7 +196,7 @@ def test_func_vs_class_KLE_FFT(): num_grid_points = ng, seed = seed) - x_t_array_class = stoc_proc.x_for_initial_time_grid() + x_t_array_class = stoc_proc.get_z() print(np.max(np.abs(x_t_array_func - x_t_array_class))) assert np.all(x_t_array_func == x_t_array_class), "stochastic_process_fft vs. StocProc Class not identical" @@ -214,7 +215,7 @@ def test_stochastic_process_KLE_interpolation(plot=False): ng_fine = ng*3 seed = 0 - sig_min = 1e-2 + sig_min = 1e-5 stoc_proc = sp.StocProc.new_instance_by_name(name = 'trapezoidal', r_tau = r_tau, @@ -229,13 +230,12 @@ def test_stochastic_process_KLE_interpolation(plot=False): ns = 6000 x_t_samples = np.empty(shape=(ns, ng_fine), dtype=np.complex) - + + print("generate samples ...") for n in range(ns): - if n % 100 == 0: - print(n, ns) stoc_proc.new_process() x_t_samples[n] = stoc_proc(finer_t) - + print("done!") ac_kle_int_conj, ac_kle_int_not_conj = sp.auto_correlation(x_t_samples) t_grid = np.linspace(0, t_max, ng_fine) @@ -293,6 +293,93 @@ def test_stochastic_process_KLE_interpolation(plot=False): assert max_diff_not_conj < 4e-2 assert max_diff_conj < 4e-2 +def test_stocproc_KLE_splineinterpolation(plot=False): + s_param = 1 + gamma_s_plus_1 = gamma(s_param+1) + # two parameter correlation function -> correlation matrix + r_tau = lambda tau : corr(tau, s_param, gamma_s_plus_1) + # time interval [0,T] + t_max = 15 + # number of subintervals + # leads to N+1 grid points + ng_fredholm = 60 + ng_kle_interp = ng_fredholm*3 + ng_fine = ng_fredholm*30 + + seed = 0 + sig_min = 1e-5 + stoc_proc = sp.StocProc_KLE(r_tau, t_max, ng_kle_interp, ng_fredholm, seed, sig_min) + + finer_t = np.linspace(0, t_max, ng_fine) + + ns = 6000 + + ac_conj = np.zeros(shape=(ng_fine, ng_fine), dtype=np.complex) + ac_not_conj = np.zeros(shape=(ng_fine, ng_fine), dtype=np.complex) + print("generate samples ...") + for n in range(ns): + stoc_proc.new_process() + x_t = stoc_proc(finer_t) + ac_conj += x_t.reshape(ng_fine, 1) * np.conj(x_t.reshape(1, ng_fine)) + ac_not_conj += x_t.reshape(ng_fine, 1) * x_t.reshape(1, ng_fine) + print("done!") + ac_conj /= ns + ac_not_conj /= ns + + t_grid = np.linspace(0, t_max, ng_fine) + ac_true = r_tau(t_grid.reshape(ng_fine, 1) - t_grid.reshape(1, ng_fine)) + + max_diff_conj = np.max(np.abs(ac_true - ac_conj)) + print("max diff : {:.2e}".format(max_diff_conj)) + + max_diff_not_conj = np.max(np.abs(ac_not_conj)) + print("max diff : {:.2e}".format(max_diff_not_conj)) + + if plot: + v_min_real = np.floor(np.min(np.real(ac_true))) + v_max_real = np.ceil(np.max(np.real(ac_true))) + + 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=3, ncols=3, figsize=(10,12)) + 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) + 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) + + ax[1,0].set_title(r"KLE $\mathrm{re}\left(\langle x(t) x^\ast(s) \rangle\right)$") + ax[1,0].imshow(np.real(ac_conj), interpolation='none', vmin=v_min_real, vmax=v_max_real) + ax[1,1].set_title(r"KLE $\mathrm{im}\left(\langle x(t) x^\ast(s) \rangle\right)$") + ax[1,1].imshow(np.imag(ac_conj), interpolation='none', vmin=v_min_imag, vmax=v_max_imag) + + ax[2,0].set_title(r"KLE $\mathrm{re}\left(\langle x(t) x(s) \rangle\right)$") + ax[2,0].imshow(np.real(ac_not_conj), interpolation='none', vmin=v_min_real, vmax=v_max_real) + ax[2,1].set_title(r"KLE $\mathrm{im}\left(\langle x(t) x(s) \rangle\right)$") + ax[2,1].imshow(np.imag(ac_not_conj), interpolation='none', vmin=v_min_imag, vmax=v_max_imag) + + ax[0,2].set_title(r"FFT log rel diff") + im02 = ax[0,2].imshow(np.log10(np.abs(ac_conj - ac_true) / np.abs(ac_true)), interpolation='none') + divider02 = make_axes_locatable(ax[0,2]) + cax02 = divider02.append_axes("right", size="10%", pad=0.05) + cbar02 = plt.colorbar(im02, cax=cax02) + + ax[1,2].set_title(r"FFT rel diff") + im12 = ax[1,2].imshow(np.abs(ac_conj - ac_true) / np.abs(ac_true), interpolation='none') + divider12 = make_axes_locatable(ax[1,2]) + cax12 = divider12.append_axes("right", size="10%", pad=0.05) + cbar12 = plt.colorbar(im12, cax=cax12) + + ax[2,2].set_title(r"FFT abs diff") + im22 = ax[2,2].imshow(np.abs(ac_conj - ac_true), interpolation='none') + divider22 = make_axes_locatable(ax[2,2]) + cax22 = divider22.append_axes("right", size="10%", pad=0.05) + cbar22 = plt.colorbar(im22, cax=cax22) + + plt.show() + + assert max_diff_not_conj < 4e-2 + assert max_diff_conj < 4e-2 def test_stochastic_process_FFT_interpolation(plot=False): s_param = 1 @@ -320,18 +407,13 @@ def test_stochastic_process_FFT_interpolation(plot=False): ac_conj = np.zeros(shape=(ng_fine, ng_fine), dtype=np.complex) ac_not_conj = np.zeros(shape=(ng_fine, ng_fine), dtype=np.complex) - - + print("generate samples ...") for n in range(ns): - if (n % 1000) == 0: - print(n, ns) stoc_proc.new_process() x_t = stoc_proc(finer_t) ac_conj += x_t.reshape(ng_fine, 1) * np.conj(x_t.reshape(1, ng_fine)) ac_not_conj += x_t.reshape(ng_fine, 1) * x_t.reshape(1, ng_fine) - - - + print("done!") ac_conj /= ns ac_not_conj /= ns @@ -479,7 +561,7 @@ def test_auto_grid_points(): # time interval [0,T] t_max = 15 ng_interpolation = 1000 - tol = 1e-16 + tol = 1e-8 ng = sp.auto_grid_points(r_tau, t_max, ng_interpolation, tol, sig_min=0) print(ng) @@ -506,7 +588,7 @@ def test_chache(): for t_i in t.keys(): for i in range(t[t_i]): total += 1 - stocproc.x(t_i) + stocproc(t_i) ci = stocproc.get_cache_info() assert ci.hits == total - misses @@ -538,25 +620,23 @@ def test_dump_load(): assert np.all(x_t == x_t_2) - - def show_auto_grid_points_result(): s_param = 1 gamma_s_plus_1 = gamma(s_param+1) # two parameter correlation function -> correlation matrix r_tau = lambda tau : corr(tau, s_param, gamma_s_plus_1) # time interval [0,T] - t_max = 3 + t_max = 15 ng_interpolation = 1000 - tol = 1e-10 + tol = 1e-8 seed = None sig_min = 0 t_large = np.linspace(0, t_max, ng_interpolation) name = 'mid_point' - name = 'trapezoidal' - name = 'gauss_legendre' +# name = 'trapezoidal' +# name = 'gauss_legendre' ng = sp.auto_grid_points(r_tau, t_max, ng_interpolation, tol, name=name, sig_min=sig_min) @@ -570,25 +650,22 @@ def show_auto_grid_points_result(): diff_max = sp._max_error(r_t_s, r_t_s_exact) plt.plot(t_large, diff) - plt.axhline(y=diff_max) + plt.plot(t_large, diff_max) plt.yscale('log') plt.grid() plt.show() - -def test_fft_func_vs_class(): - pass - - + if __name__ == "__main__": # test_stochastic_process_KLE_correlation_function(plot=False) # test_stochastic_process_FFT_correlation_function(plot=False) # test_func_vs_class_KLE_FFT() - test_stochastic_process_KLE_interpolation(plot=True) - test_stochastic_process_FFT_interpolation(plot=True) +# test_stochastic_process_KLE_interpolation(plot=False) +# test_stocproc_KLE_splineinterpolation(plot=False) +# test_stochastic_process_FFT_interpolation(plot=False) # test_stocProc_eigenfunction_extraction() # test_orthonomality() # test_auto_grid_points() -# show_auto_grid_points_result() + show_auto_grid_points_result() # test_chache() # test_dump_load() pass \ No newline at end of file