diff --git a/python/richard_hops/bcf.py b/python/richard_hops/bcf.py index 095c29b..b364262 100644 --- a/python/richard_hops/bcf.py +++ b/python/richard_hops/bcf.py @@ -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,102 +46,108 @@ 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 - - normed: if set True, scaled eta, such that int_0^infty J(w) = 1 - - mp: if set True, use mpmath for evaluation - - this class implements set_state and get_state and is suitable for binfootprint - - make sure, that s, eta, and gamma are regular float objects + represents an super / sub Ohmic spectral density function + + normed: if set True, scaled eta, such that int_0^infty J(w) = 1 + + mp: if set True, use mpmath for evaluation + + this class implements set_state and get_state and is suitable for binfootprint + + make sure, that s, eta, and gamma are regular float objects """ - + # the "args" attributes are use for get_state and set_state - self.s_arg = s - self.gamma_arg = gamma + self.s_arg = s + self.gamma_arg = gamma self.normed_arg = normed - self.mp_arg = mp + self.mp_arg = mp if normed: self.eta_arg = None else: self.eta_arg = eta - + # non non math version if not mp: self._exp = np.exp - self._gamma_func = gamma_func - self.s = s + self._gamma_func = gamma_func + self.s = s self.gamma = gamma self.eta = eta - + # mpmath version else: self._exp = mpmath.exp - self._gamma_func = mpmath.gamma - self.s = mpmath.mpf(s) + self._gamma_func = mpmath.gamma + self.s = mpmath.mpf(s) self.gamma = mpmath.mpf(gamma) self.eta = mpmath.mpf(eta) - + if normed: self.eta = self.eta / self.integral() self.c = self.eta * self._gamma_func(self.s) - + 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 + 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) + 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 def __setstate__(self, state): self.__init__(*state) + class MakeItReal(object): def __init__(self, func): self.func = func @@ -151,33 +159,31 @@ class MakeItReal(object): return self.func - - class OBCF(object): r"""ohmic bath correlation functions (BCF) - - consider ohmic spectral density - J(w) = eta w^s exp(-w / gamma) - - general BCF - alpha(tau) = 1/pi int_0^infty d w J(w) exp(-i w tau) - - -> ohmic BCF - 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) + + consider ohmic spectral density + J(w) = eta w^s exp(-w / gamma) + + general BCF + alpha(tau) = 1/pi int_0^infty d w J(w) exp(-i w tau) + + -> ohmic BCF + 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) - + if normed: self.eta = None 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(): # warnings.filterwarnings('error') @@ -192,39 +198,50 @@ 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 - + def __setstate__(self, state): eta = state[1] if eta is None: self.__init__(*state, normed=True) else: 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)" + s = "\nbcf with J(w) = eta w^s exp(-w / gamma)" s += "\ns :{}".format(self.s) s += "\neta :{}".format(self.eta) s += "\ngamma:{}".format(self.gamma) - + 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) @@ -321,76 +369,96 @@ class BCF_aprx(object): assert self._n == omega.shape self.g = g self.omega = omega - + def __call__(self, tau): r"""return alpha(tau) = sum_i=1^n g_i exp(-w_i tau) for arbitrary shape of tau""" try: 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 - + def get_SD(self, t_max, n): """ - t0, t1, t2 ... ,tn=t_max - - al(t0), al(t1), al(t2), ... al(t n-1) + al(tn), al(t n+1) = al(t n-1)^ast, ... al(t 2n-1) = al(-1)^ast - - idx: 0 ... n-1 n .. -1 - - -> dt = t_max / n-1 - + t0, t1, t2 ... ,tn=t_max + + al(t0), al(t1), al(t2), ... al(t n-1) + al(tn), al(t n+1) = al(t n-1)^ast, ... al(t 2n-1) = al(-1)^ast + + idx: 0 ... n-1 n .. -1 + + -> dt = t_max / n-1 + """ 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]) ) ) - - N = 2*(n-1) - t_, dt_ = np.linspace(0, 2*t_max, N, endpoint=False, retstep = True) - assert(dt == dt_) + 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_ assert len(t_) == len(alpha_t) - + fft_alpha = np.fft.ifft(alpha_t) * dt * N - - dw = 2*np.pi / N / dt - - freq = np.hstack( (np.arange(N//2), -np.arange(1, N//2+1)[::-1]) ) * dw + + dw = 2 * np.pi / N / dt + + 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)) - - def n(self): - return self._n[0] + return ( + self.g + * 2 + * self.omega + / ((om - 1j * self.omega) * (om + 1j * self.omega)) + ) - + def n(self): + return self._n[0] diff --git a/python/richard_hops/hierarchyData.py b/python/richard_hops/hierarchyData.py index e0bfa55..99499ce 100644 --- a/python/richard_hops/hierarchyData.py +++ b/python/richard_hops/hierarchyData.py @@ -17,18 +17,20 @@ 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" def projector_psi_t(psi_t, normed=False): """ - assume shape (len(t), dim_H_sys) + assume shape (len(t), dim_H_sys) """ psi_t_col = np.expand_dims(psi_t, axis=2) psi_t_row = np.expand_dims(psi_t, axis=1) @@ -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 + 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, - g_scale=None, - sample_method='random', - seed=0, - nonlinear=True, - normalized=False, - terminator=False, - result_type=RESULT_TYPE_ZEROTH_ORDER_ONLY, - accum_only=None, - rand_skip=None): + __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", + seed=0, + nonlinear=True, + normalized=False, + terminator=False, + result_type=RESULT_TYPE_ZEROTH_ORDER_ONLY, + accum_only=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 + 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), - 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), - dtype=np.complex128, - maxshape=(None, size_t, size_sys, size_sys), - chunks=(1, size_t, size_sys, 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( + "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), + ) + 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), - dtype=np.complex128, - maxshape=(None, size_t, size_aux_state), - chunks=(1, 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), + ) 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), - dtype=np.complex128, - maxshape=(None, size_t, num_bcf_terms), - chunks=(1, 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), + ) - 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] - if c > self.size_sys: # the linear HI has no stoc_proc data - self.stoc_proc[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 :] 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), + 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, + size_y=key.Eta.get_num_y() + ) - 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) 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['/'], - shallow=False, - expand_soft=False, - expand_external=False, - expand_refs=False, - without_attrs=False) - for k,v in f_old.attrs.items(): + 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(): 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) - - - - diff --git a/python/richard_hops/hierarchyLib.py b/python/richard_hops/hierarchyLib.py index bdcbbe2..20017f4 100644 --- a/python/richard_hops/hierarchyLib.py +++ b/python/richard_hops/hierarchyLib.py @@ -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 diff --git a/python/richard_hops/run_Fulu.py b/python/richard_hops/run_Fulu.py index 579e4dd..af0d5a1 100644 --- a/python/richard_hops/run_Fulu.py +++ b/python/richard_hops/run_Fulu.py @@ -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), - intgr_tol=1e-3, - intpl_tol=1e-3, - seed=None, - negative_frequencies=False, - scale=bcf_scale) +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, +) # 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), - intgr_tol=1e-3, - intpl_tol=1e-3, - seed=None, - negative_frequencies=False, - scale=bcf_scale) +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, +) # 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() diff --git a/python/richard_hops/run_example.py b/python/richard_hops/run_example.py index 2501d1e..76bb377 100644 --- a/python/richard_hops/run_example.py +++ b/python/richard_hops/run_example.py @@ -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), - intgr_tol=1e-2, - intpl_tol=1e-2, - seed=None, - negative_frequencies=False, - scale=coupling_strength) +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, +) # 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()