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

some_fixes

Browse source
This commit is contained in:
Valentin Boettcher 2021-10-18 18:40:50 +02:00
parent d762cc8b29
commit 46d46bbb11
5 changed files with 796 additions and 473 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

212
python/richard_hops/bcf.py
View file

@ -1,21 +1,24 @@
import numpy as np
from scipy.special import gamma as gamma_func
#from scipy.special import zeta as zeta_func
# from scipy.special import zeta as zeta_func
from scipy.integrate import quad
from collections import namedtuple
from functools import partial
#import mpmath
# import mpmath
import warnings
import traceback
import mpmath
#import mypath
#import convenientMath as cm
#import persistentdata as pd
#import binfootprint as bf
#import fcSpline
# import mypath
# import convenientMath as cm
# import persistentdata as pd
# import binfootprint as bf
# import fcSpline
import logging
log = logging.getLogger(__name__)
#############################################################
@ -24,19 +27,18 @@ log = logging.getLogger(__name__)
#
#############################################################
def zeta_func(x, q):
return np.complex128(mpmath.zeta(x, q))
obcf_param_type = namedtuple(typename = 'obcf_param_type',
field_names = ['s',
'eta',
'gamma'])
obcf_param_type = namedtuple(
typename="obcf_param_type", field_names=["s", "eta", "gamma"]
)
obcf_param_KAST_ANKERHOLD_type = namedtuple(typename = 'obcf_param_KAST_ANKERHOLD_type',
field_names = ['s',
'alpha',
'omega_c'])
obcf_param_KAST_ANKERHOLD_type = namedtuple(
typename="obcf_param_KAST_ANKERHOLD_type", field_names=["s", "alpha", "omega_c"]
)
#############################################################
#
@ -44,11 +46,13 @@ obcf_param_KAST_ANKERHOLD_type = namedtuple(typename = 'obcf_param_KAST_ANKER
#
#############################################################
class OSD(object):
r"""ohmic spectral density
J(w) = eta w^s exp(-w / gamma)
"""
def __init__(self, s, eta, gamma, normed=False, mp=False):
"""
represents an super / sub Ohmic spectral density function
@ -95,44 +99,47 @@ class OSD(object):
def __call__(self, w):
if isinstance(w, np.ndarray):
res = np.empty_like(w)
idx_l0 = np.where(w<0)
idx_l0 = np.where(w < 0)
res[idx_l0] = 0
idx_ge0 = np.where(w >= 0)
res[idx_ge0] = self.eta * w[idx_ge0]**self.s * self._exp(-w[idx_ge0]/self.gamma)
res[idx_ge0] = (
self.eta * w[idx_ge0] ** self.s * self._exp(-w[idx_ge0] / self.gamma)
)
return res
else:
if w < 0:
return 0
else:
return self.eta * w**self.s * self._exp(-w/self.gamma)
return self.eta * w ** self.s * self._exp(-w / self.gamma)
def maximum_at(self):
return self.s*self.gamma
return self.s * self.gamma
def maximum_val(self):
return self(self.maximum_at())
def integral(self):
return self.eta * self.gamma**(self.s+1) * self._gamma_func(self.s+1)
return self.eta * self.gamma ** (self.s + 1) * self._gamma_func(self.s + 1)
def reorganization_energy(self):
"""
int_0^infity J(w) / w
"""
return self.eta * self.gamma**self.s * self._gamma_func(self.s)
return self.eta * self.gamma ** self.s * self._gamma_func(self.s)
def shifted_bath_influence(self, t):
"""
int_0^infity J(w) / w exp(-i w t) = eta [ gamma/(1+i gamma t) ]**s * Gamma(s)
"""
return self.c * (self.gamma/(1+1j*self.gamma*t))**self.s
return self.c * (self.gamma / (1 + 1j * self.gamma * t)) ** self.s
def __str__(self):
return ("J(w) = eta w^s exp(-w / w_c)\n"+
"eta = {}\n".format(self.eta)+
"s = {}\n".format(self.s)+
"w_c = {}\n".format(self.gamma))
return (
"J(w) = eta w^s exp(-w / w_c)\n"
+ "eta = {}\n".format(self.eta)
+ "s = {}\n".format(self.s)
+ "w_c = {}\n".format(self.gamma)
)
def __getstate__(self):
return self.s_arg, self.eta_arg, self.gamma_arg, self.normed_arg, self.mp_arg
@ -140,6 +147,7 @@ class OSD(object):
def __setstate__(self, state):
self.__init__(*state)
class MakeItReal(object):
def __init__(self, func):
self.func = func
@ -151,8 +159,6 @@ class MakeItReal(object):
return self.func
class OBCF(object):
r"""ohmic bath correlation functions (BCF)
@ -166,6 +172,7 @@ class OBCF(object):
alpha(tau) = eta/pi gamma^(s+1) int_0^infty d x x^s exp(-x) exp(-i x gamma t)
= eta/pi gamma^(s+1) (1 + i gamma t)^(-s-1) gamma_func(s+1)
"""
def __init__(self, s, eta, gamma, normed=False):
self.s = float(s)
self.gamma = float(gamma)
@ -175,8 +182,7 @@ class OBCF(object):
self._c1 = 1 / np.pi
else:
self.eta = eta
self._c1 = self.eta * gamma_func(s + 1) * gamma**(s+1) / np.pi
self._c1 = self.eta * gamma_func(s + 1) * gamma ** (s + 1) / np.pi
def __call__(self, tau):
# with warnings.catch_warnings():
@ -192,10 +198,16 @@ class OBCF(object):
# print('a1', a1)
# print('a2', a2)
# raise
return self._c1 * (1 + 1j*self.gamma * tau)**(-(self.s+1))
return self._c1 * (1 + 1j * self.gamma * tau) ** (-(self.s + 1))
def div1(self, tau):
return - self._c1 * (self.s+1) * (1 + 1j*self.gamma*tau)**(-(self.s+2)) * 1j * self.gamma
return (
-self._c1
* (self.s + 1)
* (1 + 1j * self.gamma * tau) ** (-(self.s + 2))
* 1j
* self.gamma
)
def __getstate__(self):
return self.s, self.eta, self.gamma
@ -208,7 +220,11 @@ class OBCF(object):
self.__init__(*state, normed=False)
def __eq__(self, other):
return (self.s == other.s) and (self.eta == other.eta) and (self.gamma == other.gamma)
return (
(self.s == other.s)
and (self.eta == other.eta)
and (self.gamma == other.gamma)
)
def __repr__(self):
s = "\nbcf with J(w) = eta w^s exp(-w / gamma)"
@ -218,13 +234,14 @@ class OBCF(object):
return s
class OBCF_FT(object):
def __init__(self, s, eta, gamma, beta):
self.s_p1 = s + 1
self.eta = eta
self.gamma = gamma
self.beta = beta
self._c = self.eta / self.beta **(self.s_p1) * gamma_func(self.s_p1) / np.pi
self._c = self.eta / self.beta ** (self.s_p1) * gamma_func(self.s_p1) / np.pi
self._beta_gamma = self.beta * self.gamma
def __call__(self, tau):
@ -233,53 +250,83 @@ class OBCF_FT(object):
res_flat = res.flat
tau_flat = tau.flat
for i, ti in enumerate(tau_flat):
zf = zeta_func(self.s_p1, (1 + self._beta_gamma + 1j * self.gamma * ti) / self._beta_gamma)
res_flat[i] = self._c * ( (self._beta_gamma / (1+1j*self.gamma*ti))**self.s_p1 + zf + np.conj(zf))
zf = zeta_func(
self.s_p1,
(1 + self._beta_gamma + 1j * self.gamma * ti) / self._beta_gamma,
)
res_flat[i] = self._c * (
(self._beta_gamma / (1 + 1j * self.gamma * ti)) ** self.s_p1
+ zf
+ np.conj(zf)
)
return res
else:
zf = zeta_func(self.s_p1, (1 + self._beta_gamma + 1j * self.gamma * tau) / self._beta_gamma)
return self._c * ((self._beta_gamma / (1 + 1j * self.gamma * tau)) ** self.s_p1 + zf + np.conj(zf))
zf = zeta_func(
self.s_p1,
(1 + self._beta_gamma + 1j * self.gamma * tau) / self._beta_gamma,
)
return self._c * (
(self._beta_gamma / (1 + 1j * self.gamma * tau)) ** self.s_p1
+ zf
+ np.conj(zf)
)
def __bfkey__(self):
return self.s_p1, self.eta, self.gamma, self.beta
class OBCF_FT_mp(object):
def __init__(self, s, eta, gamma, beta):
self.s_p1 = mpmath.mpf(s) + 1
self.eta = mpmath.mpf(eta)
self.gamma = mpmath.mpf(gamma)
self.beta = mpmath.mpf(beta)
self._c = self.eta / self.beta ** (self.s_p1) * mpmath.gamma(self.s_p1) / mpmath.pi
self._c = (
self.eta / self.beta ** (self.s_p1) * mpmath.gamma(self.s_p1) / mpmath.pi
)
self._beta_gamma = self.beta * self.gamma
def __call__(self, tau):
zf = mpmath.zeta(self.s_p1, (1 + self._beta_gamma + 1j * self.gamma * tau) / self._beta_gamma)
return self._c * ((self._beta_gamma / (1 + 1j * self.gamma * tau)) ** self.s_p1 + zf + mpmath.conj(zf))
def __call__(self, tau):
zf = mpmath.zeta(
self.s_p1, (1 + self._beta_gamma + 1j * self.gamma * tau) / self._beta_gamma
)
return self._c * (
(self._beta_gamma / (1 + 1j * self.gamma * tau)) ** self.s_p1
+ zf
+ mpmath.conj(zf)
)
class OFTCorr(object):
def __init__(self, s, eta, gamma, beta):
self.s_p1 = s+1
self.s_p1 = s + 1
self.eta = eta
self.gamma = gamma
self.beta = beta
self._c = self.eta / self.beta **(self.s_p1) * gamma_func(self.s_p1) / np.pi
self._c = self.eta / self.beta ** (self.s_p1) * gamma_func(self.s_p1) / np.pi
self._beta_gamma = self.beta * self.gamma
def __call__(self, tau):
if isinstance(tau, np.ndarray):
res = np.empty(shape = tau.shape, dtype=np.complex128)
res = np.empty(shape=tau.shape, dtype=np.complex128)
res_flat = res.flat
tau_flat = tau.flat
for i, ti in enumerate(tau_flat):
res_flat[i] = self._c * zeta_func(self.s_p1, (1 + self._beta_gamma + 1j*self.gamma*ti) / self._beta_gamma)
res_flat[i] = self._c * zeta_func(
self.s_p1,
(1 + self._beta_gamma + 1j * self.gamma * ti) / self._beta_gamma,
)
return res
else:
return self._c * zeta_func(self.s_p1, (1 + self._beta_gamma + 1j*self.gamma*tau) / self._beta_gamma)
return self._c * zeta_func(
self.s_p1,
(1 + self._beta_gamma + 1j * self.gamma * tau) / self._beta_gamma,
)
def __bfkey__(self):
return self.s_p1, self.eta, self.gamma, self.beta
def BOSE_distribution_single(x):
"""calc (exp(x)-1)^-1"""
if x < 0:
@ -287,7 +334,8 @@ def BOSE_distribution_single(x):
elif x > 40:
return np.exp(-x)
else:
return 1/np.expm1(x)
return 1 / np.expm1(x)
def BOSE_distribution(x):
try:
@ -302,18 +350,18 @@ class OFTDens(object):
self.beta = beta
def __call__(self, w):
return BOSE_distribution(self.beta*w) * self.osd(w)
return BOSE_distribution(self.beta * w) * self.osd(w)
def __bfkey__(self):
return self.osd, self.beta
class BCF_aprx(object):
r"""approximation of the bath correlation function using exponentials
alpha(tau) = sum_i=1^n g_i exp(-omega_i tau)
"""
def __init__(self, g, omega):
g = np.asarray(g)
omega = np.asarray(omega)
@ -328,20 +376,22 @@ class BCF_aprx(object):
s_tau = tau.shape
dims_tau = len(s_tau)
tau_ = tau.reshape( (1,) + s_tau )
g_ = self.g.reshape( self._n + dims_tau*(1,) )
omega_ = self.omega.reshape( self._n + dims_tau*(1,) )
tau_ = tau.reshape((1,) + s_tau)
g_ = self.g.reshape(self._n + dims_tau * (1,))
omega_ = self.omega.reshape(self._n + dims_tau * (1,))
res = np.empty(shape=s_tau, dtype=np.complex128)
idx_pos = np.where(tau >= 0)
res[idx_pos] = np.sum(g_ * np.exp(-omega_*tau_[(0,)+idx_pos]), axis = 0)
res[idx_pos] = np.sum(g_ * np.exp(-omega_ * tau_[(0,) + idx_pos]), axis=0)
idx_neg = np.where(tau < 0)
res[idx_neg] = np.conj(np.sum(g_ * np.exp(-omega_*np.abs(tau_[(0,)+idx_neg])), axis = 0))
res[idx_neg] = np.conj(
np.sum(g_ * np.exp(-omega_ * np.abs(tau_[(0,) + idx_neg])), axis=0)
)
except Exception as e:
if tau >= 0:
res = np.sum(self.g*np.exp(-self.omega*tau))
res = np.sum(self.g * np.exp(-self.omega * tau))
else:
res = np.sum(self.g*np.exp(self.omega*tau)).conj()
res = np.sum(self.g * np.exp(self.omega * tau)).conj()
return res
@ -358,39 +408,57 @@ class BCF_aprx(object):
"""
t, dt = np.linspace(0, t_max, n, retstep=True)
alpha_t = self.__call__(t)
alpha_t = np.hstack( (alpha_t[:-1], np.conj(alpha_t[:0:-1]) ) )
alpha_t = np.hstack((alpha_t[:-1], np.conj(alpha_t[:0:-1])))
N = 2*(n-1)
t_, dt_ = np.linspace(0, 2*t_max, N, endpoint=False, retstep = True)
assert(dt == dt_)
N = 2 * (n - 1)
t_, dt_ = np.linspace(0, 2 * t_max, N, endpoint=False, retstep=True)
assert dt == dt_
assert len(t_) == len(alpha_t)
fft_alpha = np.fft.ifft(alpha_t) * dt * N
dw = 2*np.pi / N / dt
dw = 2 * np.pi / N / dt
freq = np.hstack( (np.arange(N//2), -np.arange(1, N//2+1)[::-1]) ) * dw
freq = np.hstack((np.arange(N // 2), -np.arange(1, N // 2 + 1)[::-1])) * dw
return 0.5 * np.fft.ifftshift(np.real(fft_alpha)), np.fft.ifftshift(freq)
def get_SD_analyt(self, w):
sd = np.asarray([np.sum((self.g.real * self.omega.real - self.g.imag*(wi-self.omega.imag)) / (self.omega.real**2 + (wi - self.omega.imag)**2)) for wi in w])
sd = np.asarray(
[
np.sum(
(
self.g.real * self.omega.real
- self.g.imag * (wi - self.omega.imag)
)
/ (self.omega.real ** 2 + (wi - self.omega.imag) ** 2)
)
for wi in w
]
)
return sd
def ft(self, om):
if isinstance(om, np.ndarray):
res = np.empty(shape = om.shape, dtype=np.complex128)
res = np.empty(shape=om.shape, dtype=np.complex128)
res_flat_view = res.flat
om_flat = om.flatten().reshape(-1, 1)
res_flat_view[:] = np.sum(self.g * 2 * self.omega / ((om_flat - 1j * self.omega) * (om_flat + 1j * self.omega)), axis=1)
res_flat_view[:] = np.sum(
self.g
* 2
* self.omega
/ ((om_flat - 1j * self.omega) * (om_flat + 1j * self.omega)),
axis=1,
)
return res
else:
return self.g * 2 * self.omega / ((om - 1j * self.omega) * (om + 1j * self.omega))
return (
self.g
* 2
* self.omega
/ ((om - 1j * self.omega) * (om + 1j * self.omega))
)
def n(self):
return self._n[0]

594
python/richard_hops/hierarchyData.py
View file

@ -17,11 +17,13 @@ import binfootprint as bf
log = logging.getLogger(__name__)
HIMetaKey_type = namedtuple('HIMetaKey_type', ['HiP', 'IntP', 'SysP', 'Eta', 'EtaTherm'])
HIMetaKey_type = namedtuple(
"HIMetaKey_type", ["HiP", "IntP", "SysP", "Eta", "EtaTherm"]
)
RESULT_TYPE_ZEROTH_ORDER_ONLY = 'ZEROTH_ORDER_ONLY'
RESULT_TYPE_ZEROTH_ORDER_AND_ETA_LAMBDA = 'ZEROTH_ORDER_AND_ETA_LAMBDA'
RESULT_TYPE_ALL = 'ALL'
RESULT_TYPE_ZEROTH_ORDER_ONLY = "ZEROTH_ORDER_ONLY"
RESULT_TYPE_ZEROTH_ORDER_AND_ETA_LAMBDA = "ZEROTH_ORDER_AND_ETA_LAMBDA"
RESULT_TYPE_ALL = "ALL"
CHAR_SET = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
@ -37,25 +39,40 @@ def projector_psi_t(psi_t, normed=False):
N = np.sum(np.conj(psi_t) * psi_t, axis=1).reshape(psi_t.shape[0], 1, 1)
return (psi_t_col * np.conj(psi_t_row)) / N
else:
return (psi_t_col * np.conj(psi_t_row))
return psi_t_col * np.conj(psi_t_row)
class HiP(object):
"""
a purely readable (non binary) data object
"""
__slots__ = ['k_max', 'g_scale', 'sample_method', 'seed',
'nonlinear', 'normalized', 'terminator', 'result_type', 'accum_only', 'rand_skip']
def __init__(self, k_max,
__slots__ = [
"k_max",
"g_scale",
"sample_method",
"seed",
"nonlinear",
"normalized",
"terminator",
"result_type",
"accum_only",
"rand_skip",
]
def __init__(
self,
k_max,
g_scale=None,
sample_method='random',
sample_method="random",
seed=0,
nonlinear=True,
normalized=False,
terminator=False,
result_type=RESULT_TYPE_ZEROTH_ORDER_ONLY,
accum_only=None,
rand_skip=None):
rand_skip=None,
):
self.k_max = k_max
self.g_scale = g_scale
self.sample_method = sample_method
@ -68,51 +85,98 @@ class HiP(object):
self.rand_skip = rand_skip
if accum_only is None:
if self.rand_skip is not None:
raise ValueError("if accum_only is 'None' (not set) rand_skip must also be 'None'")
raise ValueError(
"if accum_only is 'None' (not set) rand_skip must also be 'None'"
)
def __bfkey__(self):
if self.accum_only is None:
if self.rand_skip is not None:
raise ValueError("if accum_only is 'None' (not set) rand_skip must also be 'None'")
return [self.k_max, self.g_scale, self.sample_method,
self.seed, self.nonlinear, self.normalized, self.terminator, self.result_type]
raise ValueError(
"if accum_only is 'None' (not set) rand_skip must also be 'None'"
)
return [
self.k_max,
self.g_scale,
self.sample_method,
self.seed,
self.nonlinear,
self.normalized,
self.terminator,
self.result_type,
]
if self.rand_skip is None:
return [self.k_max, self.g_scale, self.sample_method,
self.seed, self.nonlinear, self.normalized, self.terminator, self.result_type,
self.accum_only]
return [self.k_max, self.g_scale, self.sample_method,
self.seed, self.nonlinear, self.normalized, self.terminator, self.result_type,
self.accum_only, self.rand_skip]
return [
self.k_max,
self.g_scale,
self.sample_method,
self.seed,
self.nonlinear,
self.normalized,
self.terminator,
self.result_type,
self.accum_only,
]
return [
self.k_max,
self.g_scale,
self.sample_method,
self.seed,
self.nonlinear,
self.normalized,
self.terminator,
self.result_type,
self.accum_only,
self.rand_skip,
]
def __repr__(self):
return ("k_max : {}\n".format(self.k_max)+
"g_scale : {}\n".format(self.g_scale) +
"sample_method: {}\n".format(self.sample_method) +
"seed : {}\n".format(self.seed) +
"nonlinear : {}\n".format(self.nonlinear) +
"normalized : {}\n".format(self.normalized) +
"terminator : {}\n".format(self.terminator) +
"result_type : {}\n".format(self.result_type) +
"accum_only : {}\n".format(self.accum_only) +
"rand_skip : {}\n".format(self.rand_skip) +
"")
return (
"k_max : {}\n".format(self.k_max)
+ "g_scale : {}\n".format(self.g_scale)
+ "sample_method: {}\n".format(self.sample_method)
+ "seed : {}\n".format(self.seed)
+ "nonlinear : {}\n".format(self.nonlinear)
+ "normalized : {}\n".format(self.normalized)
+ "terminator : {}\n".format(self.terminator)
+ "result_type : {}\n".format(self.result_type)
+ "accum_only : {}\n".format(self.accum_only)
+ "rand_skip : {}\n".format(self.rand_skip)
+ ""
)
class IntP(object):
"""
a purely readable (non binary) data object
"""
__slots__ = ['t_max', 't_steps', 'integrator_name', 'atol',
'rtol', 'order', 'nsteps', 'method', 't_steps_skip']
def __init__(self, t_max, t_steps,
integrator_name = 'zvode',
atol = 1e-8,
rtol = 1e-8,
order = 5,
nsteps = 5000,
method = 'bdf',
t_steps_skip = 1):
__slots__ = [
"t_max",
"t_steps",
"integrator_name",
"atol",
"rtol",
"order",
"nsteps",
"method",
"t_steps_skip",
]
def __init__(
self,
t_max,
t_steps,
integrator_name="zvode",
atol=1e-8,
rtol=1e-8,
order=5,
nsteps=5000,
method="bdf",
t_steps_skip=1,
):
self.t_max = t_max
self.t_steps = t_steps
self.integrator_name = integrator_name
@ -125,31 +189,74 @@ class IntP(object):
def __bfkey__(self):
if self.t_steps_skip == 1:
return [self.t_max, self.t_steps, self.integrator_name, self.atol,
self.rtol, self.order, self.nsteps, self.method]
return [
self.t_max,
self.t_steps,
self.integrator_name,
self.atol,
self.rtol,
self.order,
self.nsteps,
self.method,
]
else:
return [self.t_max, self.t_steps, self.integrator_name, self.atol,
self.rtol, self.order, self.nsteps, self.method, self.t_steps_skip]
return [
self.t_max,
self.t_steps,
self.integrator_name,
self.atol,
self.rtol,
self.order,
self.nsteps,
self.method,
self.t_steps_skip,
]
def __repr__(self):
return ("t_max : {}\n".format(self.t_max)+
"t_steps : {}\n".format(self.t_steps) +
"integrator_name: {}\n".format(self.integrator_name) +
"atol : {}\n".format(self.atol) +
"rtol : {}\n".format(self.rtol) +
"order : {}\n".format(self.order) +
"nsteps : {}\n".format(self.nsteps) +
"method : {}\n".format(self.method) +
"t_steps_skip : {}\n".format(self.t_steps_skip))
return (
"t_max : {}\n".format(self.t_max)
+ "t_steps : {}\n".format(self.t_steps)
+ "integrator_name: {}\n".format(self.integrator_name)
+ "atol : {}\n".format(self.atol)
+ "rtol : {}\n".format(self.rtol)
+ "order : {}\n".format(self.order)
+ "nsteps : {}\n".format(self.nsteps)
+ "method : {}\n".format(self.method)
+ "t_steps_skip : {}\n".format(self.t_steps_skip)
)
class SysP(object):
__slots__ = ['H_sys', 'L', 'psi0', 'g', 'w', 'H_dynamic', 'bcf_scale', 'gw_hash',
'len_gw', 'gw_info', 'T', 'T_method']
__slots__ = [
"H_sys",
"L",
"psi0",
"g",
"w",
"H_dynamic",
"bcf_scale",
"gw_hash",
"len_gw",
"gw_info",
"T",
"T_method",
]
def __init__(self, H_sys, L, psi0, g, w, H_dynamic, bcf_scale, gw_hash, len_gw,
gw_info = None, # these are only info fields
T = None, # which do not enter the key
T_method = None): # as T is included either in g/w or in EtaTherm
def __init__(
self,
H_sys,
L,
psi0,
g,
w,
H_dynamic,
bcf_scale,
gw_hash,
len_gw,
gw_info=None, # these are only info fields
T=None, # which do not enter the key
T_method=None,
): # as T is included either in g/w or in EtaTherm
self.H_sys = H_sys
self.L = L
self.psi0 = psi0
@ -163,59 +270,85 @@ class SysP(object):
self.T = T
self.T_method = T_method
if (self.gw_hash is None) and ( (self.g is None) or (self.w is None)):
if (self.gw_hash is None) and ((self.g is None) or (self.w is None)):
raise ValueError("specify either g/w or gw_hash")
def __bfkey__(self):
return [self.H_sys, self.L, self.psi0, self.g, self.w, self.H_dynamic, self.bcf_scale, self.gw_hash]
return [
self.H_sys,
self.L,
self.psi0,
self.g,
self.w,
self.H_dynamic,
self.bcf_scale,
self.gw_hash,
]
def __repr__(self):
return ("H_sys : {}\n".format(self.H_sys)+
"L : {}\n".format(self.L) +
"psi0 : {}\n".format(self.psi0) +
"g : {}\n".format(self.g) +
"w : {}\n".format(self.w) +
"H_dynamic: {}\n".format(self.H_dynamic) +
"bcf_scale: {}\n".format(self.bcf_scale) +
"gw_hash : {}\n".format(self.gw_hash))
return (
"H_sys : {}\n".format(self.H_sys)
+ "L : {}\n".format(self.L)
+ "psi0 : {}\n".format(self.psi0)
+ "g : {}\n".format(self.g)
+ "w : {}\n".format(self.w)
+ "H_dynamic: {}\n".format(self.H_dynamic)
+ "bcf_scale: {}\n".format(self.bcf_scale)
+ "gw_hash : {}\n".format(self.gw_hash)
)
RAND_STR_ASCII_IDX_LIST = list(range(48,58)) + list(range(65,91)) + list(range(97,123))
def rand_str(l = 8):
s = ''
RAND_STR_ASCII_IDX_LIST = (
list(range(48, 58)) + list(range(65, 91)) + list(range(97, 123))
)
def rand_str(l=8):
s = ""
for i in range(l):
s += chr(random.choice(RAND_STR_ASCII_IDX_LIST))
return s
HIData_default_size_stoc_traj = 10
HIData_default_size_rho_t_accum_part = 10
def is_int_power(x, b=2):
n_float = np.log(x)/np.log(b)
n_float = np.log(x) / np.log(b)
n = int(n_float)
if (b**n == x):
if b ** n == x:
return n
else:
return None
def file_does_not_exists_or_is_empty(fname):
if not os.path.exists(fname):
return True
else:
return (os.path.getsize(fname) == 0)
return os.path.getsize(fname) == 0
class HIData(object):
def __init__(self, hdf5_name, size_sys, size_t, size_aux_state=0, num_bcf_terms=0, accum_only=False, read_only=False):
def __init__(
self,
hdf5_name,
size_sys,
size_t,
size_y,
size_aux_state=0,
num_bcf_terms=0,
accum_only=False,
read_only=False,
):
self.hdf5_name = hdf5_name
self.accum_only = accum_only
if file_does_not_exists_or_is_empty(hdf5_name):
# print("file_does_not_exists_or_is_empty {} -> call init_file".format(hdf5_name))
self.init_file(size_sys, size_t, size_aux_state, num_bcf_terms)
self.init_file(size_sys, size_t, size_y, size_aux_state, num_bcf_terms)
else:
if not read_only:
try:
@ -223,22 +356,21 @@ class HIData(object):
except Exception as e:
print("test_file_version FAILED with exception {}".format(e))
r = input("to ignore the error type y")
if r != 'y':
if r != "y":
raise
# print("read_only", read_only)
if read_only:
# print("before h5py.File(self.hdf5_name, 'r', swmr=True, libver='latest')")
self.h5File = h5py.File(self.hdf5_name, 'r', swmr=True, libver='latest')
self.h5File = h5py.File(self.hdf5_name, "r", swmr=True, libver="latest")
# print("after")
# print(self.h5File)
else:
# print("before h5py.File(self.hdf5_name, 'r+', libver='latest')")
# print("open r+", self.hdf5_name)
try:
self.h5File = h5py.File(self.hdf5_name, 'r+', libver='latest')
self.h5File = h5py.File(self.hdf5_name, "r+", libver="latest")
except OSError:
print("FAILED to open h5 file '{}'".format(self.hdf5_name))
raise
@ -247,30 +379,37 @@ class HIData(object):
# print(self.h5File)
try:
self.stoc_traj = self.h5File['/stoc_traj']
self.rho_t_accum = self.h5File['/rho_t_accum']
self.rho_t_accum_part = self.h5File['/rho_t_accum_part']
self.stoc_traj = self.h5File["/stoc_traj"]
self.rho_t_accum = self.h5File["/rho_t_accum"]
self.rho_t_accum_part = self.h5File["/rho_t_accum_part"]
self.rho_t_accum_part_tracker = self.h5File["/rho_t_accum_part_tracker"]
self.samples = self.h5File['/samples']
self.largest_idx = self.h5File['/largest_idx']
self.samples = self.h5File["/samples"]
self.largest_idx = self.h5File["/largest_idx"]
self.tracker = self.h5File["/tracker"]
self.y = self.h5File["/y"]
if size_aux_state != 0:
self.aux_states = self.h5File['/aux_states']
self.aux_states = self.h5File["/aux_states"]
else:
self.aux_states = None
if 'aux_states' in self.h5File:
raise TypeError("HIData with aux_states=0 finds h5 file with /aux_states")
if "aux_states" in self.h5File:
raise TypeError(
"HIData with aux_states=0 finds h5 file with /aux_states"
)
if num_bcf_terms != 0:
self.stoc_proc = self.h5File['/stoc_proc']
self.stoc_proc = self.h5File["/stoc_proc"]
else:
self.stoc_proc = None
if 'stoc_proc' in self.h5File:
raise TypeError("HIData init FAILED: num_bcf_terms=0 but h5 file {} has /stoc_proc".format(self.hdf5_name))
if "stoc_proc" in self.h5File:
raise TypeError(
"HIData init FAILED: num_bcf_terms=0 but h5 file {} has /stoc_proc".format(
self.hdf5_name
)
)
self.time = self.h5File['/time']
self.time_set = self.h5File['/time_set']
self.time = self.h5File["/time"]
self.time_set = self.h5File["/time_set"]
except KeyError:
print("KeyError in hdf5 file '{}'".format(hdf5_name))
raise
@ -279,13 +418,14 @@ class HIData(object):
self.size_sys = size_sys
self.size_aux_state_plus_sys = size_aux_state + size_sys
self.num_bcf_terms = num_bcf_terms
self.size_y = size_y
self._idx_cnt = len(self.tracker)
self._idx_rho_t_accum_part_tracker_cnt = len(self.rho_t_accum_part_tracker)
def init_file(self, size_sys, size_t, size_aux_state=0, num_bcf_terms=0):
def init_file(self, size_sys, size_t, size_y, size_aux_state=0, num_bcf_terms=0):
with h5py.File(self.hdf5_name, 'w', libver='latest') as h5File:
with h5py.File(self.hdf5_name, "w", libver="latest") as h5File:
# mode 'x': Create file, fail if exists
if not self.accum_only:
@ -294,47 +434,72 @@ class HIData(object):
# need at least one stoch traj to show k_max convergence
size_stoc_traj = 1
#size_stoc_traj may be overwritten HIData_default_size_stoc_traj to account for accum_only
h5File.create_dataset('stoc_traj', (size_stoc_traj, size_t, size_sys),
# size_stoc_traj may be overwritten HIData_default_size_stoc_traj to account for accum_only
h5File.create_dataset(
"stoc_traj",
(size_stoc_traj, size_t, size_sys),
dtype=np.complex128,
maxshape=(None, size_t, size_sys),
chunks=(1, size_t, size_sys))
h5File.create_dataset('rho_t_accum', (size_t, size_sys, size_sys), dtype=np.complex128)
h5File.create_dataset('rho_t_accum_part', (HIData_default_size_rho_t_accum_part, size_t, size_sys, size_sys),
chunks=(1, size_t, size_sys),
)
h5File.create_dataset(
"y",
(size_stoc_traj, size_y),
dtype=np.complex128,
maxshape=(None, size_y),
chunks=(1, size_y),
)
h5File.create_dataset(
"rho_t_accum", (size_t, size_sys, size_sys), dtype=np.complex128
)
h5File.create_dataset(
"rho_t_accum_part",
(HIData_default_size_rho_t_accum_part, size_t, size_sys, size_sys),
dtype=np.complex128,
maxshape=(None, size_t, size_sys, size_sys),
chunks=(1, size_t, size_sys, size_sys))
chunks=(1, size_t, size_sys, size_sys),
)
h5File.create_dataset(
"rho_t_accum_part_tracker",
data=HIData_default_size_rho_t_accum_part * [False],
dtype="bool",
maxshape=(None,),
)
h5File.create_dataset("rho_t_accum_part_tracker",
data=HIData_default_size_rho_t_accum_part*[False],
dtype='bool',
maxshape=(None, ))
h5File.create_dataset("samples", (1,), dtype=np.uint32)
h5File.create_dataset("largest_idx", (1,), dtype=np.uint32)
h5File.create_dataset('samples', (1, ), dtype=np.uint32)
h5File.create_dataset('largest_idx', (1, ), dtype=np.uint32)
h5File.create_dataset("tracker",
data=HIData_default_size_stoc_traj*[False],
dtype='bool',
maxshape=(None, ))
h5File.create_dataset(
"tracker",
data=HIData_default_size_stoc_traj * [False],
dtype="bool",
maxshape=(None,),
)
if size_aux_state != 0:
# size_stoc_traj may be overwritten HIData_default_size_stoc_traj to account for accum_only
h5File.create_dataset('aux_states', (size_stoc_traj, size_t, size_aux_state),
print((size_stoc_traj, size_t, size_aux_state))
h5File.create_dataset(
"aux_states",
(size_stoc_traj, size_t, size_aux_state),
dtype=np.complex128,
maxshape=(None, size_t, size_aux_state),
chunks=(1, size_t, size_aux_state))
chunks=(1, size_t, size_aux_state),
)
if num_bcf_terms != 0:
# size_stoc_traj may be overwritten HIData_default_size_stoc_traj to account for accum_only
h5File.create_dataset('stoc_proc', (size_stoc_traj, size_t, num_bcf_terms),
h5File.create_dataset(
"stoc_proc",
(size_stoc_traj, size_t, num_bcf_terms),
dtype=np.complex128,
maxshape=(None, size_t, num_bcf_terms),
chunks=(1, size_t, num_bcf_terms))
chunks=(1, size_t, num_bcf_terms),
)
h5File.create_dataset('time', (size_t,), dtype=np.float64)
h5File.create_dataset('time_set', (1, ), dtype=np.bool)
h5File['/time_set'][0] = False
h5File.create_dataset("time", (size_t,), dtype=np.float64)
h5File.create_dataset("time_set", (1,), dtype=np.bool)
h5File["/time_set"][0] = False
def _test_file_version(self, hdf5_name):
p = get_process_accessing_file(hdf5_name)
@ -343,15 +508,19 @@ class HIData(object):
# since that other process has already changed it
return
with h5py.File(hdf5_name, libver='latest') as h5File:
with h5py.File(hdf5_name, libver="latest") as h5File:
try:
h5File.swmr_mode = True
except ValueError as e:
s = str(e)
if s.startswith("File superblock version should be at least 3"):
print("got Value Error with msg 'File superblock version should be at least 3' -> change h5 file to new version")
print(
"got Value Error with msg 'File superblock version should be at least 3' -> change h5 file to new version"
)
h5File.close()
change_file_version_to_latest(hdf5_name)
except Exception as e:
pass
def __enter__(self):
return self
@ -372,15 +541,24 @@ class HIData(object):
self.tracker.resize(size=(size,))
if not self.accum_only:
self.stoc_traj.resize(size=(size, self.size_t, self.size_sys))
self.y.resize(size=(size, self.size_y))
if self.aux_states is not None:
self.aux_states.resize(size=(size, self.size_t, self.size_aux_state_plus_sys - self.size_sys))
self.aux_states.resize(
size=(
size,
self.size_t,
self.size_aux_state_plus_sys - self.size_sys,
)
)
if self.stoc_proc is not None:
self.stoc_proc.resize(size=(size, self.size_t, self.num_bcf_terms))
# TODO: ask richard
# self.stoc_proc.resize(size=(size, self.size_t, self.num_bcf_terms))
pass
self._idx_cnt = size
def _inc_size(self, idx):
if self._idx_cnt <= idx:
new_idx_cnt = 2*self._idx_cnt
new_idx_cnt = 2 * self._idx_cnt
while new_idx_cnt <= idx:
new_idx_cnt *= 2
@ -388,12 +566,17 @@ class HIData(object):
def _resize_rho_t_accum_part(self, size):
self.rho_t_accum_part_tracker.resize(size=(size,))
self.rho_t_accum_part.resize(size=(size,self.size_t, self.size_sys, self.size_sys))
self.rho_t_accum_part.resize(
size=(size, self.size_t, self.size_sys, self.size_sys)
)
self._idx_rho_t_accum_part_tracker_cnt = size
def _inc_rho_t_accum_part_tracker_size(self, n):
if self._idx_rho_t_accum_part_tracker_cnt <= n:
new_idx_cnt = self._idx_rho_t_accum_part_tracker_cnt + HIData_default_size_rho_t_accum_part
new_idx_cnt = (
self._idx_rho_t_accum_part_tracker_cnt
+ HIData_default_size_rho_t_accum_part
)
while new_idx_cnt <= n:
new_idx_cnt += HIData_default_size_rho_t_accum_part
self._resize_rho_t_accum_part(new_idx_cnt)
@ -416,25 +599,28 @@ class HIData(object):
raise RuntimeError("can not get time, time has not been set yet.")
return self.time[:]
def new_samples(self, idx, psi_all, result_type, normed):
def new_samples(self, idx, psi_all, result_type, normed, y):
self._inc_size(idx)
if (not self.accum_only) or (idx == 0):
c = psi_all.shape[1]
self.y[idx] = y
if result_type == RESULT_TYPE_ZEROTH_ORDER_ONLY:
self.stoc_traj[idx] = psi_all
elif result_type == RESULT_TYPE_ZEROTH_ORDER_AND_ETA_LAMBDA:
self.stoc_traj[idx] = psi_all[:,:self.size_sys]
self.stoc_traj[idx] = psi_all[:, : self.size_sys]
if c > self.size_sys: # the linear HI has no stoc_proc data
self.stoc_proc[idx] = psi_all[:, self.size_sys:]
self.stoc_proc[idx] = psi_all[:, self.size_sys :]
elif result_type == RESULT_TYPE_ALL:
self.stoc_traj[idx] = psi_all[:, :self.size_sys]
self.aux_states[idx] = psi_all[:, self.size_sys:self.size_aux_state_plus_sys]
if c > self.size_aux_state_plus_sys: # the linear HI has no stoc_proc data
self.stoc_proc = psi_all[:, self.size_aux_state_plus_sys:]
self.stoc_traj[idx] = psi_all[:, : self.size_sys]
self.aux_states[idx] = psi_all[
:, self.size_sys : self.size_aux_state_plus_sys
]
if (
c > self.size_aux_state_plus_sys
): # the linear HI has no stoc_proc data
self.stoc_proc = psi_all[:, self.size_aux_state_plus_sys :]
n = is_int_power(self.samples[0]+1, b=2)
n = is_int_power(self.samples[0] + 1, b=2)
if n is not None:
self._inc_rho_t_accum_part_tracker_size(n)
@ -442,7 +628,9 @@ class HIData(object):
self.tracker[idx] = True
self.largest_idx[0] = max(self.largest_idx[0], idx)
self.rho_t_accum[:] += projector_psi_t(psi_all[:, :self.size_sys], normed=normed)
self.rho_t_accum[:] += projector_psi_t(
psi_all[:, : self.size_sys], normed=normed
)
self.samples[0] += 1
if n is not None:
@ -453,15 +641,19 @@ class HIData(object):
def get_rho_t(self, res=None):
if res is None:
res = np.empty(shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128)
res = np.empty(
shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128
)
self.rho_t_accum.read_direct(dest=res)
return res/self.get_samples()
return res / self.get_samples()
def get_rho_t_part(self, n):
if self.has_rho_t_accum_part(n):
return self.rho_t_accum_part[n]
else:
raise RuntimeError("rho_t_accum_part with index {} has not been chrunched yet".format(n))
raise RuntimeError(
"rho_t_accum_part with index {} has not been chrunched yet".format(n)
)
def get_samples(self):
return self.samples[0]
@ -480,7 +672,9 @@ class HIData(object):
if self.has_sample(idx):
return self.stoc_traj[idx]
else:
raise RuntimeError("sample with idx {} has not yet been chrunched".format(idx))
raise RuntimeError(
"sample with idx {} has not yet been chrunched".format(idx)
)
def get_sub_rho_t(self, idx_low, idx_high, normed, overwrite=False):
name = "{}_{}".format(int(idx_low), int(idx_high))
@ -489,29 +683,36 @@ class HIData(object):
if not name in self.h5File:
smp = 0
rho_t_accum = np.zeros(shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128)
rho_t_accum = np.zeros(
shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128
)
for i in range(idx_low, idx_high):
if self.has_sample(i):
smp += 1
rho_t_accum += projector_psi_t(self.stoc_traj[i], normed=normed)
rho_t = rho_t_accum / smp
h5_data = self.h5File.create_dataset(name, shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128)
h5_data = self.h5File.create_dataset(
name,
shape=(self.size_t, self.size_sys, self.size_sys),
dtype=np.complex128,
)
h5_data[:] = rho_t
h5_data.attrs['smp'] = smp
h5_data.attrs["smp"] = smp
else:
rho_t = np.empty(shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128)
self.h5File[name].read_direct(dest = rho_t)
smp = self.h5File[name].attrs['smp']
rho_t = np.empty(
shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128
)
self.h5File[name].read_direct(dest=rho_t)
smp = self.h5File[name].attrs["smp"]
return rho_t, smp
def rewrite_rho_t(self, normed):
smp = 0
rho_t_accum = np.zeros(shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128)
for i in range(self.largest_idx[0]+1):
rho_t_accum = np.zeros(
shape=(self.size_t, self.size_sys, self.size_sys), dtype=np.complex128
)
for i in range(self.largest_idx[0] + 1):
if self.has_sample(i):
smp += 1
rho_t_accum += projector_psi_t(self.stoc_traj[i], normed=normed)
@ -524,13 +725,15 @@ class HIData(object):
class HIMetaData(object):
def __init__(self, hid_name, hid_path):
self.name = hid_name
self.path = os.path.join(hid_path, '__'+self.name)
self.path = os.path.join(hid_path, "__" + self.name)
if os.path.exists(self.path):
if not os.path.isdir(self.path):
raise NotADirectoryError("the path '{}' exists but is not a directory".format(self.path))
raise NotADirectoryError(
"the path '{}' exists but is not a directory".format(self.path)
)
else:
os.mkdir(self.path)
self._fname = os.path.join(self.path, self.name+'.sqld')
self._fname = os.path.join(self.path, self.name + ".sqld")
self._l = 8
self.db = sqd.SqliteDict(filename=self._fname, autocommit=False)
@ -553,7 +756,9 @@ class HIMetaData(object):
full_name = os.path.join(self.path, fname)
try:
open(full_name, 'x').close() # open for exclusive creation, failing if the file already exists
open(
full_name, "x"
).close() # open for exclusive creation, failing if the file already exists
return fname
except FileExistsError:
pass
@ -574,7 +779,7 @@ class HIMetaData(object):
try:
hdf5_name = self.db[hashed_key][0]
except KeyError:
hdf5_name = self._new_rand_file_name(pre = self.name+'_', end='.h5')
hdf5_name = self._new_rand_file_name(pre=self.name + "_", end=".h5")
self.db[hashed_key] = (hdf5_name, key)
self.db.commit()
@ -591,37 +796,49 @@ class HIMetaData(object):
num_bcf_terms = key.SysP.len_gw
elif key.HiP.result_type == RESULT_TYPE_ALL:
num_bcf_terms = key.SysP.len_gw
size_aux_state = combLib.number_of_all_combinations_old(n = num_bcf_terms, k_max=key.HiP.k_max)
size_aux_state = combLib.number_of_all_combinations_old(
n=num_bcf_terms, k_max=key.HiP.k_max
)
if not key.HiP.nonlinear:
num_bcf_terms = 0
if hasattr(key.HiP, 'accum_only'):
if hasattr(key.HiP, "accum_only"):
accum_only = key.HiP.accum_only
else:
accum_only = False
return HIData(os.path.join(self.path, hdf5_name),
return HIData(
os.path.join(self.path, hdf5_name),
size_sys=key.SysP.H_sys.shape[0],
size_t=key.IntP.t_steps,
size_aux_state=size_aux_state,
num_bcf_terms=num_bcf_terms,
accum_only=accum_only,
read_only=read_only)
read_only=read_only,
size_y=key.Eta.get_num_y()
)
def get_process_accessing_file(fname):
cmd = 'lsof "{}"'.format(fname)
if not os.path.exists(fname):
raise FileNotFoundError(fname)
r = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, universal_newlines=True, encoding='utf8')
if r.stderr != '':
r = subprocess.run(
cmd,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
shell=True,
universal_newlines=True,
encoding="utf8",
)
if r.stderr != "":
log.info("command '{}' stderr:\n{}".format(cmd, r.stderr))
if r.returncode == 0:
# success
out = r.stdout.split('\n')
out = r.stdout.split("\n")
head = out[0].split()
idx_PID = head.index('PID')
idx_PID = head.index("PID")
pid_list = []
for l in out[1:]:
l = l.split()
@ -632,12 +849,14 @@ def get_process_accessing_file(fname):
return pid_list
else:
# failure, also happens when no process was found
if (r.stdout == ''):
log.info("lsof has non-zero return code and empty stdout -> assume not process has access to file")
if r.stdout == "":
log.info(
"lsof has non-zero return code and empty stdout -> assume not process has access to file"
)
return []
def get_rand_file_name(l=12, must_not_exists = True):
def get_rand_file_name(l=12, must_not_exists=True):
n = len(CHAR_SET)
while True:
fname = ""
@ -651,28 +870,31 @@ def change_file_version_to_latest(h5fname):
pid_list = get_process_accessing_file(h5fname)
if len(pid_list) > 0:
raise RuntimeError("can not change file version! the following processes have access to the file: {}".format(pid_list))
raise RuntimeError(
"can not change file version! the following processes have access to the file: {}".format(
pid_list
)
)
rand_fname = get_rand_file_name()
with h5py.File(rand_fname, 'w', libver='latest') as f_new:
with h5py.File(h5fname, 'r') as f_old:
with h5py.File(rand_fname, "w", libver="latest") as f_new:
with h5py.File(h5fname, "r") as f_old:
for i in f_old:
f_old.copy(i, f_new['/'],
f_old.copy(
i,
f_new["/"],
shallow=False,
expand_soft=False,
expand_external=False,
expand_refs=False,
without_attrs=False)
for k,v in f_old.attrs.items():
without_attrs=False,
)
for k, v in f_old.attrs.items():
f_new.attrs[k] = v
print("updated h5 file {} to latest version".format(os.path.abspath(h5fname)))
shutil.move(h5fname, h5fname+rand_fname+'.old')
shutil.move(h5fname, h5fname + rand_fname + ".old")
shutil.move(rand_fname, h5fname)
os.remove(h5fname+rand_fname+'.old')
os.remove(h5fname + rand_fname + ".old")
assert not os.path.exists(rand_fname)

7
python/richard_hops/hierarchyLib.py
View file

@ -813,7 +813,7 @@ def _integrate_hierarchy(arg, const_arg, c, m):
del stoc_temp_z
# t0, t1, N, f, args, x0, integrator, verbose, c, **kwargs
return ode_wrapper.integrate_cplx(c=c, args=args_dgl, **kwargs)
return ode_wrapper.integrate_cplx(c=c, args=args_dgl, **kwargs), z
class HI(object):
@ -1222,7 +1222,8 @@ class HI(object):
log.info("idx {} found in HiData, calculation SKIPPED".format(i))
if do_calculation:
t, psi_all, e_and_trb = _integrate_hierarchy(arg, const_arg, c_int, m_int)
(t, psi_all, e_and_trb), z = _integrate_hierarchy(arg, const_arg, c_int, m_int)
print(z)
if e_and_trb is not None:
print(e_and_trb[0])
print(e_and_trb[1])
@ -1230,7 +1231,7 @@ class HI(object):
if i == 0:
data.set_time(t, force=True)
data.new_samples(idx=i, psi_all=psi_all,
result_type=self.HiP.result_type, normed=self._normed_average)
result_type=self.HiP.result_type, normed=self._normed_average, y=z)
else:
skipped += 1

113
python/richard_hops/run_Fulu.py
View file

@ -14,100 +14,108 @@ numExpFit = 5
s = 1
wc = 0.531
with open('good_fit_data_abs_brute_force', "rb") as f:
with open("good_fit_data_abs_brute_force", "rb") as f:
good_fit_data_abs = pickle.load(f)
_, g_tilde, w_tilde = good_fit_data_abs[(numExpFit, s)]
g = 1/np.pi * gamma_func(s+1) * wc**(s+1) * np.asarray(g_tilde)
g = 1 / np.pi * gamma_func(s + 1) * wc ** (s + 1) * np.asarray(g_tilde)
w = wc * np.asarray(w_tilde)
bcf_scale = np.pi / 2 * alpha * wc**(1-s)
bcf_scale = np.pi / 2 * alpha * wc ** (1 - s)
# setup hierarchy parameters, here use a simple triangular hierarchy up to level 5
# the comments show the default values
# NOTE! the normalized version of HOPS is still buggy !!!
HierachyParam = hierarchyData.HiP(k_max = 5,
#g_scale=None,
#sample_method='random',
#seed=0,
#nonlinear=True,
#normalized=False,
#terminator=False,
#result_type= hierarchyData.RESULT_TYPE_ZEROTH_ORDER_ONLY,
#accum_only=None,
#rand_skip=None
)
HierachyParam = hierarchyData.HiP(
k_max=5,
# g_scale=None,
# sample_method='random',
# seed=0,
# nonlinear=True,
# normalized=False,
# terminator=False,
# result_type= hierarchyData.RESULT_TYPE_ZEROTH_ORDER_ONLY,
# accum_only=None,
# rand_skip=None
)
# setup the integration parameters (see scipy ode)
IntegrationParam = hierarchyData.IntP(t_max = t_max, t_steps = 500,
#integrator_name='zvode',
#atol=1e-8,
#rtol=1e-8,
#order=5,
#nsteps=5000,
#method='bdf',
#t_steps_skip=1
)
IntegrationParam = hierarchyData.IntP(
t_max=t_max,
t_steps=500,
# integrator_name='zvode',
# atol=1e-8,
# rtol=1e-8,
# order=5,
# nsteps=5000,
# method='bdf',
# t_steps_skip=1
)
# setup the system parameters
# here a simple Spin Boson system, single qubit coupled to a sub-ohmic environment
# how to get the g, w is a different story
SystemParam = hierarchyData.SysP(H_sys = np.asarray([[0,1],[1,0]]), # SIGMA X
L = np.asarray([[-1,0],[0,1]]), # SIGME Z
psi0 = np.asarray([0,1]), # excited qubit
g = np.asarray(g),
w = np.asarray(w),
H_dynamic = [],
bcf_scale = bcf_scale, # some coupling strength (scaling of the fit parameters 'g_i')
gw_hash = None, # this is used to load g,w from some database
len_gw = None)
SystemParam = hierarchyData.SysP(
H_sys=np.asarray([[0, 1], [1, 0]]), # SIGMA X
L=np.asarray([[-1, 0], [0, 1]]), # SIGME Z
psi0=np.asarray([0, 1]), # excited qubit
g=np.asarray(g),
w=np.asarray(w),
H_dynamic=[],
bcf_scale=bcf_scale, # some coupling strength (scaling of the fit parameters 'g_i')
gw_hash=None, # this is used to load g,w from some database
len_gw=None,
)
# setup the stochastic process for the subohmic SD
# for each sample, the StocProc_FFT class is seeded differently which results in a different realization.
# for the HOPS solver, this looks like a simple callable z(t) returning complex values
Eta = StocProc_FFT(spectral_density = bcf.OSD(s=s, eta=1, gamma=wc),
t_max = t_max,
alpha = bcf.OBCF(s=s, eta=1, gamma=wc),
Eta = StocProc_FFT(
spectral_density=bcf.OSD(s=s, eta=1, gamma=wc),
t_max=t_max,
alpha=bcf.OBCF(s=s, eta=1, gamma=wc),
intgr_tol=1e-3,
intpl_tol=1e-3,
seed=None,
negative_frequencies=False,
scale=bcf_scale)
scale=bcf_scale,
)
# thermal noise
EtaTherm = StocProc_FFT(spectral_density = bcf.OFTDens(s=s, eta=1, gamma=wc, beta=1/T),
t_max = t_max,
alpha = bcf.OFTCorr(s=s, eta=1, gamma=wc, beta=1/T),
EtaTherm = StocProc_FFT(
spectral_density=bcf.OFTDens(s=s, eta=1, gamma=wc, beta=1 / T),
t_max=t_max,
alpha=bcf.OFTCorr(s=s, eta=1, gamma=wc, beta=1 / T),
intgr_tol=1e-3,
intpl_tol=1e-3,
seed=None,
negative_frequencies=False,
scale=bcf_scale)
scale=bcf_scale,
)
# this guy must be digestable by the 'binfootprint' class, generating a unique binary representation of hi_key
# this allows to generate a hash value of the binary to access the produced data from some hdf5 data (see hierarchyData.py)
# NOTE! this requires that the ENTIRE sub-structure of 'hi_key' is digestable by 'binfootprint', therefore the special classes
# representing the spectral density and bath correlation function defined in bcf.py
# I developed 'binfootprit' since pickle and co. are not guaranteed to lead to unique binary data,
hi_key = hierarchyData.HIMetaKey_type(HiP = HierachyParam,
IntP = IntegrationParam,
SysP = SystemParam,
Eta = Eta,
EtaTherm = EtaTherm)
hi_key = hierarchyData.HIMetaKey_type(
HiP=HierachyParam,
IntP=IntegrationParam,
SysP=SystemParam,
Eta=Eta,
EtaTherm=EtaTherm,
)
# here the hierarchy structure it set up ...
myHierarchy = hierarchyLib.HI(hi_key,
number_of_samples=100,
desc="run a test case")
myHierarchy = hierarchyLib.HI(hi_key, number_of_samples=100, desc="run a test case")
# ... which allows to perform the integration, the data is stored in a file called 'myHierarchy.data'
# running the script again only crunches samples not done before!
myHierarchy.integrate_simple(data_name='myHierarchy.data')
myHierarchy.integrate_simple(data_name="myHierarchy.data")
# to access the data the 'hi_key' is used to find the data in the hdf5 file
with hierarchyData.HIMetaData(hid_name='myHierarchy.data', hid_path='.') as metaData:
with hierarchyData.HIMetaData(hid_name="myHierarchy.data", hid_path=".") as metaData:
with metaData.get_HIData(hi_key, read_only=True) as data:
smp = data.get_samples()
print("{} samples found in database".format(smp))
@ -117,7 +125,8 @@ with hierarchyData.HIMetaData(hid_name='myHierarchy.data', hid_path='.') as meta
# to show some dynamics
import matplotlib.pyplot as plt
plt.plot(t, rho_t[:,1,1].real, label="average ({})".format(smp))
plt.plot(t, rho_t[:, 1, 1].real, label="average ({})".format(smp))
plt.ylabel("rho_(0,0)")
plt.xlabel("time")
plt.legend()

137
python/richard_hops/run_example.py
View file

@ -4,6 +4,11 @@ import numpy as np
from stocproc.stocproc import StocProc_FFT
import bcf
import matplotlib
matplotlib.use("TkCairo", force=True)
t_max = 25
coupling_strength = 0.3
@ -11,64 +16,79 @@ coupling_strength = 0.3
# setup hierarchy parameters, here use a simple triangular hierarchy up to level 5
# the comments show the default values
# NOTE! the normalized version of HOPS is still buggy !!!
HierachyParam = hierarchyData.HiP(k_max = 5,
#g_scale=None,
#sample_method='random',
#seed=0,
#nonlinear=True,
#normalized=False,
#terminator=False,
#result_type= hierarchyData.RESULT_TYPE_ZEROTH_ORDER_ONLY,
#accum_only=None,
#rand_skip=None
)
HierachyParam = hierarchyData.HiP(
k_max=2,
# g_scale=None,
# sample_method='random',
# seed=0,
# nonlinear=True,
# normalized=False,
# terminator=False,
result_type=hierarchyData.RESULT_TYPE_ALL,
# accum_only=None,
# rand_skip=None
)
# setup the integration parameters (see scipy ode)
IntegrationParam = hierarchyData.IntP(t_max = t_max, t_steps = 500,
#integrator_name='zvode',
#atol=1e-8,
#rtol=1e-8,
#order=5,
#nsteps=5000,
#method='bdf',
#t_steps_skip=1
)
IntegrationParam = hierarchyData.IntP(
t_max=t_max,
t_steps=500,
# integrator_name='zvode',
# atol=1e-8,
# rtol=1e-8,
# order=5,
# nsteps=5000,
# method='bdf',
# t_steps_skip=1
)
# setup the system parameters
# here a simple Spin Boson system, single qubit coupled to a sub-ohmic environment
# how to get the g, w is a different story
SystemParam = hierarchyData.SysP(H_sys = np.asarray([[0,1],[1,0]]), # SIGMA X
L = np.asarray([[-1,0],[0,1]]), # SIGME Z
psi0 = np.asarray([0,1]), # excited qubit
g = np.asarray([-0.00783725-0.03065741j, # fit for sub-ohmic SD / BCF
-0.00959114-0.12684515j,
0.12327218-0.34116637j,
0.71838859-0.3019772j,
0.52298016+0.66270003j,
-0.20912288+0.13755251j]),
w = np.asarray([0.08064116+7.08122348e-03j, # fit for sub-ohmic SD / BCF
0.40868336+4.73041920e-02j,
1.18403114+2.25925466e-01j,
2.81218754+9.72139059e-01j,
5.63558761+3.41859384e+00j,
9.49687769+9.75760546e+00j]),
H_dynamic = [],
bcf_scale = coupling_strength, # some coupling strength (scaling of the fit parameters 'g_i')
gw_hash = None, # this is used to load g,w from some database
len_gw = None)
SystemParam = hierarchyData.SysP(
H_sys=np.asarray([[0, 1], [1, 0]]), # SIGMA X
L=np.asarray([[-1, 0], [0, 1]]), # SIGME Z
psi0=np.asarray([0, 1]), # excited qubit
g=np.asarray(
[
-0.00783725 - 0.03065741j, # fit for sub-ohmic SD / BCF
-0.00959114 - 0.12684515j,
0.12327218 - 0.34116637j,
0.71838859 - 0.3019772j,
0.52298016 + 0.66270003j,
-0.20912288 + 0.13755251j,
]
),
w=np.asarray(
[
0.08064116 + 7.08122348e-03j, # fit for sub-ohmic SD / BCF
0.40868336 + 4.73041920e-02j,
1.18403114 + 2.25925466e-01j,
2.81218754 + 9.72139059e-01j,
5.63558761 + 3.41859384e00j,
9.49687769 + 9.75760546e00j,
]
),
H_dynamic=[],
bcf_scale=coupling_strength, # some coupling strength (scaling of the fit parameters 'g_i')
gw_hash=None, # this is used to load g,w from some database
len_gw=None,
)
# setup the stochastic process for the subohmic SD
# for each sample, the StocProc_FFT class is seeded differently which results in a different realization.
# for the HOPS solver, this looks like a simple callable z(t) returning complex values
Eta = StocProc_FFT(spectral_density = bcf.OSD(s=0.25, eta=1, gamma=3),
t_max = t_max,
alpha = bcf.OBCF(s=0.25, eta=1, gamma=3),
Eta = StocProc_FFT(
spectral_density=bcf.OSD(s=0.25, eta=1, gamma=3),
t_max=t_max,
alpha=bcf.OBCF(s=0.25, eta=1, gamma=3),
intgr_tol=1e-2,
intpl_tol=1e-2,
seed=None,
negative_frequencies=False,
scale=coupling_strength)
scale=coupling_strength,
)
# no thermal noise -> zero temperature, however, this works ion the same fashion as above
EtaTherm = None
@ -77,36 +97,39 @@ EtaTherm = None
# NOTE! this requires that the ENTIRE sub-structure of 'hi_key' is digestable by 'binfootprint', therefore the special classes
# representing the spectral density and bath correlation function defined in bcf.py
# I developed 'binfootprit' since pickle and co. are not guaranteed to lead to unique binary data,
hi_key = hierarchyData.HIMetaKey_type(HiP = HierachyParam,
IntP = IntegrationParam,
SysP = SystemParam,
Eta = Eta,
EtaTherm = EtaTherm)
hi_key = hierarchyData.HIMetaKey_type(
HiP=HierachyParam,
IntP=IntegrationParam,
SysP=SystemParam,
Eta=Eta,
EtaTherm=EtaTherm,
)
# here the hierarchy structure it set up ...
myHierarchy = hierarchyLib.HI(hi_key,
number_of_samples=100,
desc="run a test case")
myHierarchy = hierarchyLib.HI(hi_key, number_of_samples=100, desc="run a test case")
# ... which allows to perform the integration, the data is stored in a file called 'myHierarchy.data'
# running the script again only crunches samples not done before!
myHierarchy.integrate_simple(data_name='myHierarchy.data')
myHierarchy.integrate_simple(data_name="myHierarchy.data")
# to access the data the 'hi_key' is used to find the data in the hdf5 file
with hierarchyData.HIMetaData(hid_name='myHierarchy.data', hid_path='.') as metaData:
with hierarchyData.HIMetaData(hid_name="myHierarchy.data", hid_path=".") as metaData:
with metaData.get_HIData(hi_key, read_only=True) as data:
smp = data.get_samples()
print("{} samples found in database".format(smp))
t = data.get_time()
rho_t = data.get_rho_t()
psi_zt_0 = data.get_stoc_traj(idx = 0)
psi_zt_0 = data.get_stoc_traj(idx=0)
# to show some dynamics
import matplotlib.pyplot as plt
plt.plot(t, rho_t[:,0,0].real, label="average ({})".format(smp))
plt.plot(t, np.abs(psi_zt_0[:,0]), label="single psi_z_t")
plt.plot(t, np.sqrt(np.sum(np.abs(psi_zt_0)**2, axis=(1,))), label="|psi_z_t|", color='0.4')
plt.plot(t, rho_t[:, 0, 0].real, label="average ({})".format(smp))
plt.plot(t, np.abs(psi_zt_0[:, 0]), label="single psi_z_t")
plt.plot(
t, np.sqrt(np.sum(np.abs(psi_zt_0) ** 2, axis=(1,))), label="|psi_z_t|", color="0.4"
)
plt.ylabel("rho_(0,0)")
plt.xlabel("time")
plt.legend()