hiro/stocproc
1
0
Fork 0
mirror of https://github.com/vale981/stocproc synced 2025-03-06 02:01:41 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

removed meta data

Browse source
This commit is contained in:
Richard Hartmann 2014-09-16 13:10:21 +02:00
parent 35c3ea432d
commit daa54994a5
1 changed files with 27 additions and 134 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

161
stocproc.py
View file

@ -18,8 +18,7 @@
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA 02110-1301, USA. # MA 02110-1301, USA.
""" """
Stochastic Process Module **Stochastic Process Module**
=========================
This module contains various methods to generate stochastic processes for a This module contains various methods to generate stochastic processes for a
given correlation function. There are two different kinds of generators. The one kind given correlation function. There are two different kinds of generators. The one kind
@ -57,18 +56,12 @@ solutions of the time discrete version.
.. todo:: implement convenient classes with fixed weights .. todo:: implement convenient classes with fixed weights
""" """
from __future__ import division
import numpy as np import numpy as np
import time as tm import time as tm
# _pp_ import my_pypreprocessor as pp
# _pp_ pypreprocessor = pp.Preprocessor(__file__, verbose = False, always_parse = False)
# _pp_ # pypreprocessor.define('TIME_IT')
# _pp_ pypreprocessor.define('VERBOSE')
# _pp_ # pypreprocessor.define('TESTING')
# _pp_ pypreprocessor.process()
class StocProc(object): class StocProc(object):
r"""Simulate Stochastic Process using Karhunen-Loève expansion r"""Simulate Stochastic Process using Karhunen-Loève expansion
@ -180,12 +173,6 @@ class StocProc(object):
same weights :math:`w_i` and time grid points :math:`s_i`. same weights :math:`w_i` and time grid points :math:`s_i`.
""" """
tmp = self._Y * self._sqrt_eig_val.reshape(self._num_ev,1) tmp = self._Y * self._sqrt_eig_val.reshape(self._num_ev,1)
#exclude
# _pp_ print "tmp = rand * sqrt(eig_val): shape", tmp.shape
# _pp_ print tmp
# _pp_ print "eig_vec: shape", self._eig_vec.shape
# _pp_ print self._eig_vec
#endexclude
return np.tensordot(tmp, self._eig_vec, axes=([0],[1])).flatten() return np.tensordot(tmp, self._eig_vec, axes=([0],[1])).flatten()
def x(self, t): def x(self, t):
@ -330,9 +317,7 @@ def auto_grid_points(r_tau, t_max, ng_interpolation, tol = 1e-8, err_method = _m
seed = None seed = None
sig_min = 0 sig_min = 0
t_large = np.linspace(0, t_max, ng_interpolation) t_large = np.linspace(0, t_max, ng_interpolation)
#ifdef VERBOSE print("start auto_grid_points, determine ng ...")
print "start auto_grid_points, determine ng ..."
#endif
#exponential increase to get below error threshold #exponential increase to get below error threshold
while err > tol: while err > tol:
@ -345,19 +330,17 @@ def auto_grid_points(r_tau, t_max, ng_interpolation, tol = 1e-8, err_method = _m
err = np.max(err_method(r_t_s, r_t_s_exact)) err = np.max(err_method(r_t_s, r_t_s_exact))
#ifdef VERBOSE print(" ng {} -> err {:.2e}".format(ng, err))
print " ng {} -> err {:.2e}".format(ng, err)
#endif
ng_low = ng // 2 ng_low = ng // 2
ng_high = ng ng_high = ng
while (ng_high - ng_low) > 1: while (ng_high - ng_low) > 1:
print "#"*40 print("#"*40)
print " ng_l", ng_low print(" ng_l", ng_low)
print " ng_h", ng_high print(" ng_h", ng_high)
ng = (ng_low + ng_high) // 2 ng = (ng_low + ng_high) // 2
print " ng", ng print(" ng", ng)
t, w = get_trapezoidal_weights_times(t_max, ng) t, w = get_trapezoidal_weights_times(t_max, ng)
stoc_proc = StocProc(r_tau, t, w, seed, sig_min) stoc_proc = StocProc(r_tau, t, w, seed, sig_min)
@ -366,16 +349,14 @@ def auto_grid_points(r_tau, t_max, ng_interpolation, tol = 1e-8, err_method = _m
err = np.max(err_method(r_t_s, r_t_s_exact)) err = np.max(err_method(r_t_s, r_t_s_exact))
#ifdef VERBOSE print(" ng {} -> err {:.2e}".format(ng, err))
print " ng {} -> err {:.2e}".format(ng, err)
#endif
if err > tol: if err > tol:
print " err > tol!" print(" err > tol!")
print " ng_l -> ", ng print(" ng_l -> ", ng)
ng_low = ng ng_low = ng
else: else:
print " err <= tol!" print(" err <= tol!")
print " ng_h -> ", ng print(" ng_h -> ", ng)
ng_high = ng ng_high = ng
@ -416,17 +397,13 @@ def solve_hom_fredholm(r, w, eig_val_min):
d = np.diag(np.sqrt(w)) d = np.diag(np.sqrt(w))
d_inverse = np.diag(1/np.sqrt(w)) d_inverse = np.diag(1/np.sqrt(w))
r = np.dot(d, np.dot(r, d)) r = np.dot(d, np.dot(r, d))
#ifdef VERBOSE print("solve eigenvalue equation ...")
print "solve eigenvalue equation ..."
#endif
eig_val, eig_vec = np.linalg.eigh(r) eig_val, eig_vec = np.linalg.eigh(r)
# use only eigenvalues larger than sig_min**2 # use only eigenvalues larger than sig_min**2
large_eig_val_idx = np.where(eig_val >= eig_val_min)[0] large_eig_val_idx = np.where(eig_val >= eig_val_min)[0]
num_of_functions = len(large_eig_val_idx) num_of_functions = len(large_eig_val_idx)
#ifdef VERBOSE print("use {} / {} eigenfunctions".format(num_of_functions, len(w)))
print "use {} / {} eigenfunctions".format(num_of_functions, len(w))
#endif
eig_val = eig_val[large_eig_val_idx] eig_val = eig_val[large_eig_val_idx]
eig_vec = eig_vec[:, large_eig_val_idx] eig_vec = eig_vec[:, large_eig_val_idx]
@ -489,10 +466,8 @@ def stochastic_process_kle(r_tau, t, w, num_samples, seed = None, sig_min = 1e-4
stochastic process. stochastic process.
""" """
#ifdef VERBOSE print("__ stochastic_process __")
print "__ stochastic_process __" print("pre calculations ...")
print "pre calculations ..."
#endif
if seed != None: if seed != None:
np.random.seed(seed) np.random.seed(seed)
@ -503,9 +478,6 @@ def stochastic_process_kle(r_tau, t, w, num_samples, seed = None, sig_min = 1e-4
# r_tau(t-s) -> integral/sum over s -> s must be row in EV equation # r_tau(t-s) -> integral/sum over s -> s must be row in EV equation
r = r_tau(t_col-t_row) r = r_tau(t_col-t_row)
#ifdef TIME_IT
# _pp_ t_0 = tm.clock()
#endif
# solve discrete Fredholm equation # solve discrete Fredholm equation
# eig_val = lambda # eig_val = lambda
@ -515,61 +487,14 @@ def stochastic_process_kle(r_tau, t, w, num_samples, seed = None, sig_min = 1e-4
# generate samples # generate samples
sig = np.sqrt(eig_val).reshape(1, num_of_functions) # variance of the random quantities of the Karhunen-Loève expansion sig = np.sqrt(eig_val).reshape(1, num_of_functions) # variance of the random quantities of the Karhunen-Loève expansion
#ifdef TIME_IT
# _pp_ t_1 = tm.clock()
#endif
#ifdef VERBOSE print("generate samples ...")
print "generate samples ..."
#endif
x = np.random.normal(size=(num_samples, num_of_functions)) * sig # random quantities all aligned for num_samples samples x = np.random.normal(size=(num_samples, num_of_functions)) * sig # random quantities all aligned for num_samples samples
#exclude
# _pp_ print "x = rand * sqrt(eigval): shape", x.shape
# _pp_ print x
# _pp_ print "eig_vec: shape", eig_vec.shape
# _pp_ print eig_vec
#endexclude
x_t_array = np.tensordot(x, eig_vec, axes=([1],[1])) # multiplication with the eigenfunctions == base of Karhunen-Loève expansion x_t_array = np.tensordot(x, eig_vec, axes=([1],[1])) # multiplication with the eigenfunctions == base of Karhunen-Loève expansion
#ifdef TESTING
# _pp_ print "use np.dot ..."
# _pp_ x_t_array_2 = np.dot(x, eig_vec.T)
# _pp_ assert np.all(x_t_array == x_t_array_2)
# _pp_ print "passed test: np.all(x_t_array == x_t_array_2)"
#endif
#ifdef TIME_IT
# _pp_ t_2 = tm.clock()
#endif
#ifdef TESTING
#ifdef VERBOSE
# _pp_ print "generate samples (alternative) ..."
#endif
# _pp_ assert seed != None
# _pp_ x_t_array_alt = _stochastic_process_alternative_samples(num_samples, num_of_functions, t, sig.flatten(), eig_vec, seed)
#ifdef TIME_IT
# _pp_ t_3 = tm.clock()
#endif
# _pp_ assert np.all(x_t_array == x_t_array_alt)
# _pp_ print "passed test: assert np.all(x_t_array == x_t_array_alt)"
#ifdef TIME_IT
# _pp_ print "generate samples alternative : {:.2e}s".format(t_3-t_2)
# _pp_ print "-"*50
#endif
#endif
#ifdef TIME_IT
# _pp_ print "-"*50
# _pp_ print "solve eigenvalue equation : {:.2e}s".format(t_1-t_0)
# _pp_ print "generate samples : {:.2e}s".format(t_2-t_1)
# _pp_ print "-"*50
# _pp_ print "time per sample : {:.2e}s".format((t_2-t_1) / num_samples )
# _pp_ print "overall time : {:.2e}s".format(t_2-t_0)
# _pp_ print "-"*50
#endif
#ifdef VERBOSE print("ALL DONE!\n")
print "ALL DONE!\n"
#endif
return x_t_array return x_t_array
@ -582,7 +507,7 @@ def _stochastic_process_alternative_samples(num_samples, num_of_functions, t, si
""" """
np.random.seed(seed) np.random.seed(seed)
x_t_array = np.empty(shape=(num_samples, len(t)), dtype = complex) x_t_array = np.empty(shape=(num_samples, len(t)), dtype = complex)
for i in xrange(num_samples): for i in range(num_samples):
x = np.random.normal(size=num_of_functions) * sig x = np.random.normal(size=num_of_functions) * sig
x_t_array[i,:] = np.dot(eig_vec, x.T) x_t_array[i,:] = np.dot(eig_vec, x.T)
return x_t_array return x_t_array
@ -737,13 +662,8 @@ def stochastic_process_fft(spectral_density, t_max, num_grid_points, num_samples
1D array of time grid points). Each row of the stochastic process 1D array of time grid points). Each row of the stochastic process
array contains one sample of the stochastic process. array contains one sample of the stochastic process.
""" """
#ifdef VERBOSE print("__ stochastic_process_fft __")
print "__ stochastic_process_fft __" print("pre calculations ...")
print "pre calculations ..."
#endif
#ifdef TIME_IT
# _pp_ t_0 = tm.clock()
#endif
n_dft = num_grid_points * 2 - 1 n_dft = num_grid_points * 2 - 1
delta_t = t_max / (num_grid_points-1) delta_t = t_max / (num_grid_points-1)
delta_omega = 2 * np.pi / delta_t / n_dft delta_omega = 2 * np.pi / delta_t / n_dft
@ -752,47 +672,20 @@ def stochastic_process_fft(spectral_density, t_max, num_grid_points, num_samples
omega = delta_omega*np.arange(n_dft) omega = delta_omega*np.arange(n_dft)
#reshape for multiplication with matrix xi #reshape for multiplication with matrix xi
sqrt_spectral_density = np.sqrt(spectral_density(omega)).reshape((1, n_dft)) sqrt_spectral_density = np.sqrt(spectral_density(omega)).reshape((1, n_dft))
#ifdef TIME_IT
# _pp_ t_1 = tm.clock()
#endif
if seed != None: if seed != None:
np.random.seed(seed) np.random.seed(seed)
#ifdef VERBOSE print(" omega_max : {:.2}".format(delta_omega * n_dft))
print " omega_max : {:.2}".format(delta_omega * n_dft) print(" delta_omega: {:.2}".format(delta_omega))
print " delta_omega: {:.2}".format(delta_omega) print("generate samples ...")
print "generate samples ..."
#endif
#ifdef TIME_IT
# _pp_ t_2 = tm.clock()
#endif
#random normal samples #random normal samples
xi = np.random.normal(size = (num_samples,n_dft)) xi = np.random.normal(size = (num_samples,n_dft))
#each row contain a different integrand #each row contain a different integrand
weighted_integrand = sqrt_spectral_density * xi * np.sqrt(delta_omega) weighted_integrand = sqrt_spectral_density * xi * np.sqrt(delta_omega)
#compute integral using fft routine #compute integral using fft routine
z_ast = np.fft.rfft(weighted_integrand, axis = 1) z_ast = np.fft.rfft(weighted_integrand, axis = 1)
#ifdef TIME_IT
# _pp_ t_3 = tm.clock()
#endif
#ifdef TIME_IT
# _pp_ print "-"*50
# _pp_ print "inital array setup : {:.2e}s".format(t_1-t_0)
# _pp_ print "generate samples using fft : {:.2e}s".format(t_3-t_2)
# _pp_ print "-"*50
# _pp_ print "time per sample : {:.2e}s".format( (t_3-t_2) / num_samples )
# _pp_ print "overall time : {:.2e}s".format(t_3-t_0)
# _pp_ print "-"*50
#endif
#corresponding time axis #corresponding time axis
t = np.linspace(0, t_max, num_grid_points) t = np.linspace(0, t_max, num_grid_points)
#ifdef TESTING print("ALL DONE!\n")
# _pp_ t_ = np.fft.rfftfreq(n_dft, delta_omega/2/np.pi)
# _pp_ assert (np.max(np.abs(t-t_)) < 1e-14)
# _pp_ print "passed test: assert np.all(t == t_)"
#endif
#ifdef VERBOSE
print "ALL DONE!\n"
#endif
return z_ast, t return z_ast, t
@ -820,5 +713,5 @@ def auto_correlation(x, s_0_idx = 0):
return np.mean(x * np.conj(x_s_0), axis = 0) return np.mean(x * np.conj(x_s_0), axis = 0)
if __name__ == '__main__': if __name__ == '__main__':
print 'this is the stocproc module' print('this is the stocproc module')