mirror of
https://github.com/vale981/stocproc
synced 2025-03-04 17:21:42 -05:00
added align_eig_vec feature for testing, which fixes the overall phase of the eigen vectors for KLE
This commit is contained in:
Richard Hartmann
parent
f73c034f91
commit
6d2c6a0f1b
3 changed files with 71 additions and 17 deletions
|
|
@ -133,7 +133,7 @@ class StocProc_KLE(_absStocProc):
|
|||
- 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, verbose=1, k=3):
|
||||
def __init__(self, r_tau, t_max, ng_fredholm, ng_fac=4, seed=None, sig_min=1e-5, verbose=1, 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]
|
||||
|
|
@ -154,7 +154,8 @@ class StocProc_KLE(_absStocProc):
|
|||
ng = ng_fredholm,
|
||||
sig_min = sig_min,
|
||||
seed = seed,
|
||||
verbose = verbose)
|
||||
verbose = verbose,
|
||||
align_eig_vec = align_eig_vec)
|
||||
|
||||
ng = ng_fac * (ng_fredholm - 1) + 1
|
||||
|
||||
|
|
|
|||
|
|
@ -118,7 +118,8 @@ class StocProc(object):
|
|||
sig_min = 1e-4,
|
||||
fname = None,
|
||||
cache_size = 1024,
|
||||
verbose = 1):
|
||||
verbose = 1,
|
||||
align_eig_vec = False):
|
||||
|
||||
self.verbose = verbose
|
||||
self._one_over_sqrt_2 = 1/np.sqrt(2)
|
||||
|
|
@ -144,6 +145,12 @@ class StocProc(object):
|
|||
# eig_val = lambda
|
||||
# eig_vec = u(t)
|
||||
self._eig_val, self._eig_vec = solve_hom_fredholm(r, w, sig_min**2, verbose=self.verbose)
|
||||
if align_eig_vec:
|
||||
for i in range(self._eig_vec.shape[1]):
|
||||
s = np.sum(self._eig_vec[:,i])
|
||||
phase = np.exp(1j*np.arctan2(np.real(s), np.imag(s)))
|
||||
self._eig_vec[:,i]/= phase
|
||||
|
||||
else:
|
||||
self.__load(fname)
|
||||
|
||||
|
|
@ -154,17 +161,17 @@ class StocProc(object):
|
|||
self.new_process(seed = seed)
|
||||
|
||||
@classmethod
|
||||
def new_instance_by_name(cls, name, r_tau, t_max, ng, seed, sig_min, verbose=1):
|
||||
def new_instance_by_name(cls, name, r_tau, t_max, ng, seed, sig_min, verbose=1, align_eig_vec=False):
|
||||
known_names = ['trapezoidal', 'mid_point', 'simpson', 'gauss_legendre']
|
||||
|
||||
if name == 'trapezoidal':
|
||||
ob = cls.new_instance_with_trapezoidal_weights(r_tau, t_max, ng, seed, sig_min, verbose)
|
||||
ob = cls.new_instance_with_trapezoidal_weights(r_tau, t_max, ng, seed, sig_min, verbose, align_eig_vec)
|
||||
elif name == 'mid_point':
|
||||
ob = cls.new_instance_with_mid_point_weights(r_tau, t_max, ng, seed, sig_min, verbose)
|
||||
ob = cls.new_instance_with_mid_point_weights(r_tau, t_max, ng, seed, sig_min, verbose, align_eig_vec)
|
||||
elif name == 'simpson':
|
||||
ob = cls.new_instance_with_simpson_weights(r_tau, t_max, ng, seed, sig_min, verbose)
|
||||
ob = cls.new_instance_with_simpson_weights(r_tau, t_max, ng, seed, sig_min, verbose, align_eig_vec)
|
||||
elif name == 'gauss_legendre':
|
||||
ob = cls.new_instance_with_gauss_legendre_weights(r_tau, t_max, ng, seed, sig_min, verbose)
|
||||
ob = cls.new_instance_with_gauss_legendre_weights(r_tau, t_max, ng, seed, sig_min, verbose, align_eig_vec)
|
||||
else:
|
||||
raise RuntimeError("unknown name '{}' to create StocProc instance\nknown names are {}".format(name, known_names))
|
||||
|
||||
|
|
@ -172,24 +179,24 @@ class StocProc(object):
|
|||
return ob
|
||||
|
||||
@classmethod
|
||||
def new_instance_with_trapezoidal_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1):
|
||||
def new_instance_with_trapezoidal_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1, align_eig_vec=False):
|
||||
t, w = get_trapezoidal_weights_times(t_max, ng)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose, align_eig_vec=align_eig_vec)
|
||||
|
||||
@classmethod
|
||||
def new_instance_with_simpson_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1):
|
||||
def new_instance_with_simpson_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1, align_eig_vec=False):
|
||||
t, w = get_simpson_weights_times(t_max, ng)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose, align_eig_vec=align_eig_vec)
|
||||
|
||||
@classmethod
|
||||
def new_instance_with_mid_point_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1):
|
||||
def new_instance_with_mid_point_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1, align_eig_vec=False):
|
||||
t, w = get_mid_point_weights(t_max, ng)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose, align_eig_vec=align_eig_vec)
|
||||
|
||||
@classmethod
|
||||
def new_instance_with_gauss_legendre_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1):
|
||||
def new_instance_with_gauss_legendre_weights(cls, r_tau, t_max, ng, seed=None, sig_min=0, verbose=1, align_eig_vec=False):
|
||||
t, w = gquad.gauss_nodes_weights_legendre(n=ng, low=0, high=t_max)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose)
|
||||
return cls(r_tau, t, w, seed, sig_min, verbose=verbose, align_eig_vec=align_eig_vec)
|
||||
|
||||
def __load(self, fname):
|
||||
with open(fname, 'rb') as f:
|
||||
|
|
|
|||
|
|
@ -1134,6 +1134,51 @@ def test_ac_vs_ac_from_c():
|
|||
assert np.max(np.abs(ac - ac_c)) < 1e-15
|
||||
assert np.max(np.abs(ac_prime - ac_prime_c)) < 1e-15
|
||||
|
||||
def test_align_eig_vecs():
|
||||
r_tau = lambda tau: (1 + 1j*(tau))**(-1.8)
|
||||
|
||||
t_max = 30
|
||||
ng_fredholm = 61
|
||||
ng_fac = 4
|
||||
num_grid_points = 100
|
||||
seed = 0
|
||||
verbose = 0
|
||||
sigmin = 0
|
||||
|
||||
t_fine = np.linspace(0, t_max, num_grid_points * 10 + 1)
|
||||
|
||||
scale = 2.4
|
||||
r_tau_scaled = lambda tau: scale * r_tau(tau)
|
||||
|
||||
for align_eig_vec in [False, True]:
|
||||
|
||||
# init sp class
|
||||
sp_kle_scl = sp.class_stocproc.StocProc_KLE(r_tau=r_tau_scaled,
|
||||
t_max=t_max,
|
||||
ng_fredholm=ng_fredholm,
|
||||
ng_fac=ng_fac,
|
||||
seed=seed,
|
||||
sig_min=sigmin,
|
||||
verbose=verbose,
|
||||
align_eig_vec=align_eig_vec)
|
||||
sp_kle = sp.class_stocproc.StocProc_KLE(r_tau=r_tau,
|
||||
t_max=t_max,
|
||||
ng_fredholm=ng_fredholm,
|
||||
ng_fac=ng_fac,
|
||||
seed=seed,
|
||||
sig_min=sigmin,
|
||||
verbose=verbose,
|
||||
align_eig_vec=align_eig_vec)
|
||||
|
||||
# if eig_vecs are aligned this ration should be a constant
|
||||
if align_eig_vec:
|
||||
assert(np.allclose(sp_kle_scl._A / sp_kle._A, 1/np.sqrt(scale)))
|
||||
else:
|
||||
assert not np.allclose(sp_kle_scl._A / sp_kle._A, 1 / np.sqrt(scale))
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
|
@ -1164,5 +1209,6 @@ if __name__ == "__main__":
|
|||
# test_integral_equation()
|
||||
|
||||
# show_auto_grid_points_result()
|
||||
# show_ef()
|
||||
# show_ef()
|
||||
test_align_eig_vecs()
|
||||
pass
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue