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
master-thesis/python/richard_hops/hierarchyLib.py

1365 lines
42 KiB
Python
Raw Blame History

from collections import namedtuple
from functools import wraps, partial
import logging
import numpy as np
import scipy.sparse as spe
import sobolLib as sbl
from time import time, sleep
import warnings
import combLib as cl
import hierarchyData as hid
import ode_wrapper
import progression as progress
from numpy import concatenate as np_concatenate
from numpy import vdot as np_vdot
from numpy import empty_like as np_empty_like
from numpy import asarray as np_asarray
from numpy import sum as np_sum
from numpy import real as np_real
T_DIFF_TOL = 1e-13
log = logging.getLogger(__name__)
def time_it(f):
@wraps(f)
def wrapper(*args, **kwds):
t0 = time()
res = f(*args, **kwds)
t1 = time()
print("{}: {:.2e}s".format(f.__name__, t1 - t0))
return res
return wrapper
def get_O_vector(omega, n, kmax, idxDict):
"""
accounts for the coupling via
set notation
sum_(l=1)^(k) omega_(lambda_l) psi(t)^(lambda_1,...lambda_k)
vector notation
< K | Omega > psi(t)^(K) = (k_0 Omega_0 + k_1 Omega_1 + ... k_n Omega_n) psi(t)^[k_0,k_1,... k_n]
"""
size = len(idxDict)
o_vec = np.empty(shape=(size, 1), dtype=np.complex128) # row vector
for binkey in idxDict.keys():
w = 0
k = cl.binkey_to_nparray(binkey)
w = sum(k * omega)
o_vec[idxDict[binkey], 0] = -w
return o_vec
def get_M_up(g, n, idxDict, g_scale=None):
"""
g: list of g_i, which occur in the BCF alpha(tau) = sum_i g_i exp(-omega_i tau)
for finite bath corresponds to the coupling constant
n: number of such g_i, respectively the number of exponentials in the BCF
for finite bath corresponds to the number of oscillators
kmax: hierarchy depth
idxDict: mapping of the hierarchy function signature (set of indices) to
a single index addressing the position of the HF in a vector
g_scale: account for the freedom we have to define the HF
for each g_i there can be an arbitrary complex g_scale_i
accounts for
sum_lambda g_lambda/g_scale_lambda psi(t)^(lambda_1,...lambda_k, lambda)
sum_lambda g_lambda/g_scale_lambda psi(t)^(k + e_\lambda)
we are looking at the hierarchy equation with signature ckmo (called reference)
its time derivative will be influenced by hierarchy equations
of higher order (reference coupled "to" that equation)
so in final matrix notation
d / dt psi = M psi
the index to the reference function (idx_ref) is row
where as the idx_to is the column index
"""
assert len(g) == n
if g_scale is None:
g_scale = np.ones(n)
else:
assert len(g_scale) == n
data = []
row = []
col = []
for k in idxDict:
# each kappa as np array
k_npa = cl.binkey_to_nparray(k)
# the corresponding index
idx_ref = idxDict[k]
# up coupling, from below to hierarchy level k
for i in range(n):
# increase each entry in kappa by one
k_to_npa = np.copy(k_npa)
k_to_npa[i] += 1
# the corresponding bytes as key
k_to = k_to_npa.data.tobytes()
# if the key macthes an equation considered
# in the truncated hierarchy
if k_to in idxDict:
idx_to = idxDict[k_to]
data.append(g[i] / g_scale[i])
row.append(idx_ref)
col.append(idx_to)
size = len(idxDict)
return spe.coo_matrix((data, (row, col)), shape=(size, size)).tocsr()
def get_M_down(n, idxDict, g_scale=None):
"""
n: number of such g_i (see above), respectively the number of exponentials in the BCF
for finite bath corresponds to the number of oscillators
kmax: hierarchy depth
idxDict: mapping of the hierarchy function signature (set of indices) to
a single index addressing the position of the HF in a vector
g_scale: account for the freedom we have to define the HF
for each g_i there can be an arbitrary complex g_scale_i
accounts for
sum_(l=1)^k g_scale_lambda_l psi(t)^(lambda_1,...lambda_k)/lambda
sum_lambda kappa_lambda g_scale_lambda_l psi(t)^(k - e_\lambda)
we are looking at the hierarchy equation with signature ckmo (called reference)
its time derivative will be influenced by hierarchy equations
of higher order (reference coupled "to" that equation)
so in final matrix notation
d / dt psi = M psi
the index to the reference function (idx_ref) is row
where as the idx_to is the column index
"""
if g_scale is None:
g_scale = np.ones(n)
else:
assert len(g_scale) == n
row = []
col = []
data = []
for k in idxDict:
# each kappa as np array
k_npa = cl.binkey_to_nparray(k)
# the corresponding index
idx_ref = idxDict[k]
# up coupling, from below to hierarchy level k
for i in range(n):
# increase each entry in kappa by one
if k_npa[i] == 0:
continue
k_to_npa = np.copy(k_npa)
k_to_npa[i] -= 1
# the corresponding bytes as key
k_to = k_to_npa.data.tobytes()
idx_to = idxDict[k_to]
data.append(k_npa[i] * g_scale[i])
row.append(idx_ref)
col.append(idx_to)
size = len(idxDict)
return spe.coo_matrix((data, (row, col)), shape=(size, size)).tocsr()
F_Args_type = namedtuple(
"F_Args_type", ["eta", "O_vec", "K", "B", "C", "M_down", "M_up", "spn"]
)
def x_to_res_zeroth_order_only(t_, x_, dim_sys):
return x_[:dim_sys]
def x_to_res_zeroth_order_and_eta_lambda(t_, x_, dim_sys, num_bcf_terms):
return np_concatenate((x_[:dim_sys], x_[-num_bcf_terms:]))
def f_hierarchy(t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn):
"""
general hierarchy form
d/dt Psi = O x Psi + (K(t) x Psi^T)^T
+ (B(t) x (M_up x Psi)^T)^T
+ (C(t) x (M_down x Psi)^T)^T
note that this function should not be directly called by the integration routine
because eta_det is not respected by the integration call from hierarchy lib
so the various implementations of the hierarchy set eta_det in different
manners, in order to actually calculate the derivative, and provide
the correct function needed for he integration.
see: f_lin_hierarchy, f_non_lin_hierarchy
"""
# O.shape = (num_hier, 1)
num_hier = O_vec.shape[0]
psi = V_psi.reshape((num_hier, spn.dim_sys))
# assert np_all(V_psi[:spn.dim_sys] == psi[0,:])
# this is the result vector
ddt_psi = np_empty_like(V_psi)
# create a view on the result vector ...
ddt_psi_view = ddt_psi.view()
# ... let's us changing the shape while operating on the
# same block of memory
ddt_psi_view.shape = (num_hier, spn.dim_sys)
# filling the data here ...
K_mat = K(t, psi, spn, eta_stoc, eta_det)
ddt_psi_view[:, :] = (
O_vec * psi
+ K_mat.dot(psi.T).T
+ B(t, psi, spn).dot(M_up.dot(psi).T).T
+ C(t, psi, spn).dot(M_down.dot(psi).T).T
)
# ... leads to the correct linear alignment for the result vector
return ddt_psi
def f_hierarchy_timed(t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn):
"""
general hierarchy form
d/dt Psi = O x Psi + (K(t) x Psi^T)^T
+ (B(t) x (M_up x Psi)^T)^T
+ (C(t) x (M_down x Psi)^T)^T
note that this function should not be directly called by the integration routine
because eta_det is not respected by the integration call from hierarchy lib
so the various implementations of the hierarchy set eta_det in different
manners, in order to actually calculate the derivative, and provide
the correct function needed for he integration.
see: f_lin_hierarchy, f_non_lin_hierarchy
"""
# O.shape = (num_hier, 1)
t0 = time()
num_hier = O_vec.shape[0]
psi = V_psi.reshape((num_hier, spn.dim_sys))
# this is the result vector
ddt_psi = np_empty_like(V_psi)
# create a view on the result vector ...
ddt_psi_view = ddt_psi.view()
# ... let's us changing the shape while operating on the
# same block of memory
ddt_psi_view.shape = (num_hier, spn.dim_sys)
# filling the data here ...
t_init = time() - t0
t0 = time()
K_ = K(t, psi, spn, eta_stoc, eta_det)
t_K = time() - t0
t0 = time()
B_ = B(t, psi, spn)
t_B = time() - t0
t0 = time()
C_ = C(t, psi, spn)
t_C = time() - t0
t0 = time()
K_psi = K_.dot(psi.T).T
t_K_psi = time() - t0
t0 = time()
M_up_psi = M_up.dot(psi).T
t_M_up_psi = time() - t0
t0 = time()
M_down_psi = M_down.dot(psi).T
t_M_down_psi = time() - t0
t0 = time()
B_psi = B_.dot(M_up_psi).T
t_B_psi = time() - t0
t0 = time()
C_psi = C_.dot(M_down_psi).T
t_C_psi = time() - t0
t0 = time()
O_vec_psi = O_vec * psi
t_O_psi = time() - t0
t0 = time()
psi = O_vec_psi + K_psi + B_psi + C_psi
t_sum_psi = time() - t0
total = (
t_init
+ t_K
+ t_B
+ t_C
+ t_K_psi
+ t_M_up_psi
+ t_M_down_psi
+ t_B_psi
+ t_C_psi
+ t_O_psi
+ t_sum_psi
)
print("f_hierarchy_timed total {:.2g}".format(total))
print(" init {:.2g} ({:%})".format(t_init, t_init / total))
print(" K {:.2g} ({:%})".format(t_K, t_K / total))
print(" B {:.2g} ({:%})".format(t_B, t_B / total))
print(" C {:.2g} ({:%})".format(t_C, t_C / total))
print(" K_psi {:.2g} ({:%})".format(t_K_psi, t_K_psi / total))
print(" M_up_psi {:.2g} ({:%})".format(t_M_up_psi, t_M_up_psi / total))
print(" M_down_psi {:.2g} ({:%})".format(t_M_down_psi, t_M_down_psi / total))
print(" B_psi {:.2g} ({:%})".format(t_B_psi, t_B_psi / total))
print(" C_psi {:.2g} ({:%})".format(t_C_psi, t_C_psi / total))
print(" O_psi {:.2g} ({:%})".format(t_O_psi, t_O_psi / total))
print(" sum_psi {:.2g} ({:%})".format(t_sum_psi, t_sum_psi / total))
ddt_psi_view[:, :] = psi
# ... leads to the correct linear alignment for the result vector
return ddt_psi
def f_lin_hierarchy(t, V_psi, eta_stoc, O_vec, K, B, C, M_down, M_up, spn):
"""
this wrapper implements the linear hierarchy
eta_det is not used, and simply set to None
V_psi is a vector of size (num_hierarchies_eq * dim_sys)
"""
eta_det = None
return f_hierarchy(t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn)
def f_lin_hierarchy_timed(t, V_psi, eta_stoc, O_vec, K, B, C, M_down, M_up, spn):
eta_det = None
return f_hierarchy_timed(
t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn
)
def f_non_lin_hierarchy(
t, V_psi_and_eta_lambda, eta_stoc, O_vec, K, B, C, M_down, M_up, spn
):
"""
this wrapper calculates the derivative of the composed vector
(psi, eta_lambda)
where psi hold the hierarchy functions as in the linear case
and eta_lambda the deterministic part of the process eta(t)
"""
V_psi = V_psi_and_eta_lambda[: -spn.num_bcf_terms]
eta_lambda = V_psi_and_eta_lambda[-spn.num_bcf_terms :]
eta_det = eta_det_non_lin(spn, eta_lambda)
ddt_eta_lambda = ddt_eta_lambda_non_lin(spn, V_psi, eta_lambda)
ddt_psi = f_hierarchy(
t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn
)
res = np_empty_like(V_psi_and_eta_lambda)
res[: -spn.num_bcf_terms] = ddt_psi
res[-spn.num_bcf_terms :] = ddt_eta_lambda
return res
def f_non_lin_hierarchy_timed(
t, V_psi_and_eta_lambda, eta_stoc, O_vec, K, B, C, M_down, M_up, spn
):
t0 = time()
V_psi = V_psi_and_eta_lambda[: -spn.num_bcf_terms]
eta_lambda = V_psi_and_eta_lambda[-spn.num_bcf_terms :]
t_sep = time() - t0
t0 = time()
eta_det = eta_det_non_lin(spn, eta_lambda)
t_eta = time() - t0
t0 = time()
ddt_eta_lambda = ddt_eta_lambda_non_lin(spn, V_psi, eta_lambda)
t_ddt_eta = time() - t0
t0 = time()
ddt_psi = f_hierarchy(
t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn
)
t_ddt_psi = time() - t0
t0 = time()
res = np_empty_like(V_psi_and_eta_lambda)
res[: -spn.num_bcf_terms] = ddt_psi
res[-spn.num_bcf_terms :] = ddt_eta_lambda
t_concat = time() - t0
total = t_sep + t_eta + t_ddt_eta + t_ddt_psi + t_concat
print("f_non_lin_hierarchy_timed total {:.2e}".format(total))
print(" sep {:.2e} ({:%})".format(t_sep, t_sep / total))
print(" eta_det {:.2e} ({:%})".format(t_eta, t_eta / total))
print(" ddt eta {:.2e} ({:%})".format(t_ddt_eta, t_ddt_eta / total))
print(" ddt psi {:.2e} ({:%})".format(t_ddt_psi, t_ddt_psi / total))
print(" concat {:.2e} ({:%})".format(t_concat, t_concat / total))
return res
################################################################
##
## LINAR HIERARCHY
##
################################################################
def K_lin(t, psi, spn, eta_stoc, eta_det):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
minus_i_H_t = spn.minus_i_H_sys.copy()
for hdyn in spn.H_dyn:
minus_i_H_t += -1j * hdyn(t)
return minus_i_H_t + np.conj(eta_stoc(t)) * spn.L
def B_lin(t, psi, spn):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
return spn.minus_L_dagger
def C_lin(t, psi, spn):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
return spn.L
################################################################
##
## NON-LINEAR HIERARCHY
##
################################################################
def E(spn, psi_zero):
"""
psi_zero = psi[0:dim_sys]
E(L_dagger, psi_zero) = <psi_zero|L_dagger|psi_zero> / <psi_zero|psi_zero>
"""
res = np_vdot(psi_zero, np_asarray(spn.L_dagger.dot(psi_zero))) / np_vdot(
psi_zero, psi_zero
)
return res
def K_non_lin(t, psi, spn, eta_stoc, eta_det):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
minus_i_H_t = spn.minus_i_H_sys.copy()
for hdyn in spn.H_dyn:
minus_i_H_t += -1j * hdyn(t)
eta_stoc_t = np.conj(eta_stoc(t))
return minus_i_H_t + (eta_stoc_t + eta_det) * spn.L
def B_non_lin(t, psi, spn):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
psi_zero = psi[0, :]
return spn.minus_L_dagger + E(spn, psi_zero) * spn.eye
def C_non_lin(t, psi, spn):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
return spn.L
def eta_det_non_lin(spn, eta_lambda_t):
return np_sum(spn.g_ast * eta_lambda_t)
def ddt_eta_lambda_non_lin(spn, V_psi, eta_lambda):
"""
V_psi is a vector of size (num_hierarchies_eq * dim_sys)
"""
psi_zero = V_psi[: spn.dim_sys]
return E(spn, psi_zero) - spn.omega_ast * eta_lambda
################################################################
##
## NON-LINEAR HIERARCHY NORMALIZED
##
################################################################
def E_normalized(spn, psi_zero):
"""
psi_zero = psi[0:dim_sys]
E(L_dagger, psi_zero) = <psi_zero|L_dagger|psi_zero>
"""
res = np_vdot(psi_zero, np_asarray(spn.L_dagger.dot(psi_zero)))
return res
def K_non_lin_normalized(t, psi, spn, eta_stoc, eta_det):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
eta_t_L = (np.conj(eta_stoc(t)) + eta_det) * spn.L
psi_zero = psi[0, :]
n_exps = len(spn.g_num)
psi_one = psi[1 : n_exps + 1, :]
tmp = 0.0
E_tmp = spn.minus_L_dagger + E_normalized(spn, psi_zero) * spn.eye
for i in range(n_exps):
tmp += np_real(spn.g_num[i] * np_vdot(psi_zero, E_tmp.dot(psi_one[i, :])))
return (
spn.minus_i_H_sys
+ eta_t_L
+ (tmp - np_vdot(psi_zero, ((eta_t_L + eta_t_L.getH()) / 2).dot(psi_zero)))
* spn.eye
)
# return spn.minus_i_H_sys + eta_t_L
def K_non_lin_normalized_timed(t, psi, spn, eta_stoc, eta_det):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
t0 = time()
eta_t = eta_stoc(t)
t_eval_eta = time() - t0
t0 = time
eta_t_L = (eta_t + eta_det) * spn.L
t_eta_t_L = time() - t0
t0 = time()
psi_zero = psi[0, :]
n_exps = len(spn.g_num)
psi_one = psi[1 : n_exps + 1, :]
t_init = time() - t0
t0 = time()
tmp = 0.0
E_tmp = spn.minus_L_dagger + E_normalized(spn, psi_zero) * spn.eye
for i in range(n_exps):
tmp += np_real(spn.g_num[i] * np_vdot(psi_zero, E_tmp.dot(psi_one[i, :])))
return (
spn.minus_i_H_sys
+ eta_t_L
+ (tmp - np_vdot(psi_zero, ((eta_t_L + eta_t_L.getH()) / 2).dot(psi_zero)))
* spn.eye
)
# return spn.minus_i_H_sys + eta_t_L
def B_non_lin_normalized(t, psi, spn):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
psi_zero = psi[0, :]
return spn.minus_L_dagger + E_normalized(spn, psi_zero) * spn.eye
def C_non_lin_normalized(t, psi, spn):
"""
NOTE: psi has shape (num_hierarchies_eq, dim_sys)
"""
return spn.L
class H_stoc_pot(object):
def __init__(self, sp, L):
assert isinstance(L, np.ndarray)
self.L = L
self.L_dagger = np.conj(self.L.T)
if sp is not None:
self.sp = sp
self.new_process(np.zeros(self.get_num_y()))
else:
self.sp = lambda t: 0
def new_process(self, z):
self.sp.new_process(z)
def get_num_y(self):
return self.sp.get_num_y()
def __call__(self, t):
spt = self.sp(t)
return spt * self.L_dagger + np.conj(spt) * self.L
def __getstate__(self):
return self.L, self.sp
def __setstate__(self, state):
self.L, self.sp = state
self.L_dagger = np.conj(self.L.T)
class SysParamForNumerics(object):
def __init__(self, H_sys, L, H_dynamic, g, w, g_scale, EtaTherm=None):
if hasattr(H_sys, "todense"):
self.minus_i_H_sys = np.asarray(-1j * H_sys.todense(), dtype=np.complex128)
else:
self.minus_i_H_sys = np.asarray(-1j * H_sys, dtype=np.complex128)
if hasattr(L, "todense"):
self.L = np.asarray(L.todense(), dtype=np.complex128)
else:
self.L = np.asarray(L, dtype=np.complex128)
if hasattr(H_dynamic, "__call__"):
# assume there is a single time dependent H
self.H_dyn = [H_dynamic]
else:
# assume an iterable containing time dependent H
for h in H_dynamic:
assert hasattr(h, "__call__")
self.H_dyn = list(H_dynamic)
self.L_dagger = np.conj(self.L.T)
self.minus_L_dagger = -self.L_dagger
self.dim_sys = self.L.shape[0]
self.num_bcf_terms = len(g)
self.omega_ast = np.conj(w)
self.eye = np.eye(self.L.shape[0])
if g_scale is None:
g_scale = 1
self.g_num = g / g_scale
self.g_ast = np.conj(g) # this is not scaled, occurs in eta_det_non_lin
if EtaTherm is not None:
self.H_dyn.insert(0, H_stoc_pot(sp=EtaTherm, L=self.L))
class EtaFinitebathForNumerics(object):
"""
for the hamiltonian
H(t) = H_sys + L sum_l gamma_l^ast exp(omega_l t) a_l^dagger + L^dagger sum_l gamma_l exp(omega_l t) a_l
with pure imaginary omega,
we define the stochastic process as
eta(t) = i sum_l gamma_l exp(-omega_l t) z_l
with z_l being a gaussian complex random variable with <z z^ast> = 1
with the properties
<eta(t)> = 0
<eta(t) eta(s)> = 0
<eta(t) eta^ast(s)> = sum_l |gamma_l|^2 exp(-i |omega_l| (t-s)) == sum_l g_l exp(-i |omega_l| (t-s))
"""
def __init__(self, gamma, omega, T=None):
self.z_ast = None
self._gamma = np.asarray(gamma)
self._omega = np.asarray(omega)
self._T = T
if not (np.ndim(self._gamma) == 1):
raise RuntimeError("g must be 1D")
if not (self._omega.shape == self._gamma.shape):
raise RuntimeError("omega and gamma must have the same length")
if not np.all(np.real(self._omega) == 0):
raise RuntimeError("for finite bath, omega must be purely imaginary")
if not np.all(np.imag(self._gamma) == 0):
raise RuntimeError("for finite bath, gamma should be real")
if self._T is not None:
if self._T == 0:
bar_n = 0
else:
bar_n = np.exp(-np.abs(self._omega) / self._T) / (
1 - np.exp(-np.abs(self._omega) / self._T)
)
self._gamma = np.sqrt(bar_n) * self._gamma
def __repr__(self):
return (
"EtaFinitebathForNumerics (stochastic process class)\n"
+ "gamma: {}\n".format(self._gamma)
+ "omega: {}\n".format(self._omega)
+ "T : {} (for T not None: effective gamma = sqrt(bar n(T)) gamma)".format(
self._T
)
)
def new_process(self, z):
self.z = z
def __call__(self, t):
return 1j * np.sum(self._gamma * np.exp(-self._omega * t) * self.z)
def __getstate__(self):
return (self._gamma, self._omega)
def __setstate__(self, state):
self.z_ast = None
self._gamma, self._omega = state
def get_num_y(self):
return len(self._omega)
class Eta_ZERO_ForNumerics(object):
def __init__(self):
pass
def new_process(self, z):
pass
def __call__(self, t):
return 0
def __getstate__(self):
return 0
def __setstate__(self, state):
pass
def get_num_y(self):
return 0
args_keyword_as_named_tupled = namedtuple(
"args_keyword_as_named_tupled", ["id", "sub_seed"]
)
def merge_arg_and_const_arg(arg, const_arg):
"""
prepares data from arg and const_arg such that they can be passed
to the general integration routine
arg and const_arg are both assumed to be dictionaries
the merge process must not alter arg nor const_arg
in order to be used in the jobmanager context
returns the arguments passed to the function
defining the derivative such that
args_dgl = arg['args'] + const_arg['args']
where as arg['args'] and const_arg['args'] have been assumed to be tuples
e.g.
arg['args'] = (2, 'low')
const_arg['args'] = (15, np.pi)
f will be called with
f(t, x, 2, 'low', 15, np.pi)
returns further the combined dictionary
arg + const_arg with the keyword 'args' removed
For any duplicate keys the value will be the value
from the 'arg' dictionary.
"""
# allows arg to be a namedtuple (or any other object that
# can be converted to a dict)
# the purpose is that dicts are not allowed as keys for
# the persistent data structure, where as namedtupled are
# and therefore arg as a neamedtuple may be used to identify
# the result of this calculation in the database
if hasattr(arg, "_asdict"):
arg = arg._asdict()
# same for const_arg, as a reason of consistency
if hasattr(const_arg, "_asdict"):
const_arg = const_arg._asdict()
# extract the args keyword from arg and const_arg
args_dgl = tuple()
if "args" in arg:
args_dgl += arg["args"]
if "args" in const_arg:
args_dgl += const_arg["args"]
kwargs = {}
kwargs.update(const_arg)
kwargs.update(arg)
# remove args as they have been constructed explicitly
if "args" in kwargs:
del kwargs["args"]
# remove id, when it comes from the persistentDataServer
if "id" in kwargs:
del kwargs["id"]
return args_dgl, kwargs
def _integrate_hierarchy(arg, const_arg, c, m):
"""
arg: named tuple (id, subseed)
"""
# named tupled to dict conversion moved to merge_arg_and_const_arg
# args_dgl is simple tuple
# removes also 'id'
args_dgl, kwargs = merge_arg_and_const_arg(arg, const_arg)
m.value = kwargs["N"]
np.random.seed(arg.sub_seed)
rand_skip = kwargs["rand_skip"]
np.random.rand(rand_skip)
z_num = kwargs["z_num"]
stoc_temp_z_num = kwargs["stoc_temp_z_num"]
kwargs.pop("z_num")
kwargs.pop("stoc_temp_z_num")
kwargs.pop("sub_seed")
kwargs.pop("rand_skip")
z = uni_to_gauss(np.random.rand(z_num))
# this is a stoc proc instance, we set the process according to the RV z
args_dgl[0].new_process(z)
if stoc_temp_z_num != 0:
stoc_temp_z = uni_to_gauss(np.random.rand(stoc_temp_z_num))
args_dgl[-1].H_dyn[-1].new_process(stoc_temp_z)
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), z
class HI(object):
"""
setup efficient integration of hierarchy system of ODE
the general form is such that
d/dt psi^(i_1, .. i_k) = omega^(i_1, .. i_k) psi^(i_1, .. i_k)
+ K(t) psi^(i_1, .. i_k)
+ B(t) sum_(j=1)^n G_j / g_j psi^(i_1, .. i_k, j)
+ C(t) sum_(l=1)^k g_(i_l) psi^(i_1, .. i_k) / i_l
the 'signature' (i_1, .. i_k) addresses a unique N dimensional function
where as N corresponds to the dimension of the reduced Hilbert space
if there are k indices in the signature we call psi^(i_1, .. i_k)
a HF of order k
each index can take value between 1 ... n, where n gives the number
of exponentials in the BCF (number of oscillators for a finite bath)
so for n=2 and k=3 we have the signatures
(1,1,1) (1,1,2) (1,2,2) (2,2,2)
note: the order of the indices does not matter, rather think of
a set of indices. with that, a division (such as (i_1, .. i_k) / i_l) is
also explained by taking the i_l th index out of the set.
a unique mapping of the signatures to an index is achieved as follows
() -> 0
(1) -> 1
(2) -> 2
(1,1) -> 3
(1,2) -> 4
(2,2) -> 5
(1,1,1) -> 6
(1,1,2) -> 7
(1,2,2) -> 8
(2,2,2) -> 9
...
note, that the mapping depends crucially on n (the number of different indices)
alternatively the set can be cast to a vector of dimension corresponding to the
number of different indices (number of exponential summands in the bcf). Each
entry in the vector represents one "index type" where the value in each entry
tell the occurences of this "index type"
eg. suppose we have 2 index types, namely 1 and 2
then the set (1, 2) corresponds to the vector [1, 1] because each index type occures only onece
() == [0,0]
(1,1) == [2,0]
(1,2,2) = [1,2] and so on.
to use scipy ODE we have to put all psi in one huge vector, called V_psi,
using the index mapping above.
V_psi looks like the following
[ psi^()_1, psi^()_2, ... psi^()_N,
psi^(1)_1, psi^(1)_2, ... psi^(1)_N,
psi^(2)_1, psi^(2)_2, ... psi^(2)_N,
psi^(1,1)_1, psi^(1,1)_2, ... psi^(1,1)_N, and so on ]
to to the summation for the time derivative, we reshape V_phi -> Psi such that
Psi = [ [ psi^(1)_1, psi^(1)_2, ... psi^(1)_N ],
[ psi^(2)_1, psi^(2)_2, ... psi^(2)_N ],
[ psi^(1,1)_1, psi^(1,1)_2, ... psi^(1,1)_N],
[...],
...]
so the index mapping tells in which row to find the HF with given signature
a sparse matrix is generated, such that a simple matrix multiplication
gives the index dependent self-, up-, down couplings
(coupling HF of order k to k, k+1, k-1)
this matrix includes the corresponding index dependencies
for eg. let's call M_up the matrix which handles the up coupling then M_up
has the following structure. considering the i-th row, means looking at
the coupling of the i-th HF. the cell (i,j) then tells how the time
derivative of the HF i depends on the HF j.
for example in the case of n=2 the second row (i=1) im M_up looks like that
i:0 [ ... ],
1 [0, 0, 0, G_1/g_1, G_2/g_2]
2 [ ... ],
which leads to M_up x Psi = [ [ ...],
[ G_1/g_1*psi^(1,1)_1 + G_2/g_2*psi^(1,2)_1,
G_1/g_1*psi^(1,1)_2 + G_2/g_2*psi^(1,2)_2,
...
G_1/g_1*psi^(1,1)_N + G_2/g_2*psi^(1,2)_N],
[ ...],
...]
so M_up x Psi gives a matrix, where each row holds the correct element wise
addition of HF influencing the time derivative "from above".
Depending on the time and kind of the hierarchy each row (as thought of a
vector) will additionally modified via matrix multiplication with B(t) where
B(t) has the shape (NxN). As we have the dimension of the reduced Hilbert space
aligned in row manner we have to transpose M_up x Psi.
Psi_up = B(t) x (M_up x Psi)^T
we manage to express the up coupling of the form
B(t) sum_(j=1)^n G_j / g_j psi^(i_1, .. i_k, j)
for an individual signature (i_1, .. i_k), in matrix form for all signatures
simultaneously.
The same follows for the self and down coupling.
The task is now to control the various kinds of matrices A, B, C (see notes).
above A was split in '-sum Omega + K(t)' because Omega depends on the signature
where as K(t) is the same for all HF.
summary
-------
in matrix form we find
d/dt Psi = O x Psi + K(t) x Psi^T
+ B(t) x (M_down x Psi)^T
+ C(t) x (M_up x Psi)^T
where we have for the above example
O = [O^(), O^(1), O^(2), O^(1,1), ...] (row vector)
with O^(i_1, .. i_k) = sum_(l=1)^k omega_{i_l}
and omega_i are the n 'complex frequencies' occurring in the exponentials
of the BCF
once we have set up O, M_down, M_up at init time
we can integrate the system by providing the time dependent
(and eventually state dependent) matrices K,B,C
further we distinguish between matrix and scalar components occurring in
K, B, C as the scalar components correspond to simply multiply the matrix
Psi with a scalar value
"""
def __init__(
self, hi_key, number_of_samples, desc, min_sample_index=0, fbcf_data=None
):
self.hi_key = hi_key
self.HiP, self.IntP, self.SysP, self.Eta, self.EtaTherm = hi_key
assert isinstance(self.SysP, hid.SysP)
assert isinstance(self.IntP, hid.IntP)
assert isinstance(self.HiP, hid.HiP)
self.rand_skip = 0
if self.HiP.rand_skip is not None:
self.rand_skip = self.HiP.rand_skip
self.Eta.set_scale(self.SysP.bcf_scale)
if self.EtaTherm is not None:
self.EtaTherm.set_scale(self.SysP.bcf_scale)
self.min_sample_index = min_sample_index
self.number_of_samples = number_of_samples
self.desc = desc
self._normed_average = self.HiP.nonlinear and not self.HiP.normalized
# get g/w from data base
if self.SysP.g is None or self.SysP.w is None:
fit_data = fbcf_data.get_result(gw_hash=self.SysP.gw_hash)
self._g = fit_data.g * self.SysP.bcf_scale
self._w = fit_data.w
fbcf_data.close()
else:
self._g = self.SysP.g * self.SysP.bcf_scale
self._w = self.SysP.w
self._n = len(self._g)
warnings.warn(
"sum_k_max is not implemented! DO SO BEFORE NEXT USAGE (use simplex)."
+ "HierarchyParametersType does not yet know about sum_k_max"
)
self.idxDict = cl.idx_dict(n=self._n, k_max=self.HiP.k_max)
self.O_vec = get_O_vector(
omega=self._w, n=self._n, kmax=self.HiP.k_max, idxDict=self.idxDict
)
self.M_up = get_M_up(
g=self._g, n=self._n, idxDict=self.idxDict, g_scale=self.HiP.g_scale
)
self.M_down = get_M_down(
n=self._n, idxDict=self.idxDict, g_scale=self.HiP.g_scale
)
self.spn = SysParamForNumerics(
H_sys=self.SysP.H_sys,
L=self.SysP.L,
H_dynamic=self.SysP.H_dynamic,
g=self._g,
w=self._w,
g_scale=self.HiP.g_scale,
EtaTherm=self.EtaTherm,
)
self.numEq = len(self.idxDict)
self._dim_sys = self.SysP.H_sys.shape[0]
self._dim = self.numEq * self._dim_sys
print("init Hi class, use {} equation".format(self._dim))
self.x0 = None
if not self.HiP.nonlinear:
# simple linear case
self.K = K_lin
self.B = B_lin
self.C = C_lin
self.f_hierarchy = f_lin_hierarchy
self.f_hi_timed = f_lin_hierarchy_timed
self.x0 = self.get_x0_zeroTemp()
else:
if not self.HiP.normalized:
# non linear case
self.K = K_non_lin
self.B = B_non_lin
self.C = C_non_lin
self.f_hierarchy = f_non_lin_hierarchy
self.f_hi_timed = f_non_lin_hierarchy_timed
self.x0 = np.concatenate(
(
self.get_x0_zeroTemp(),
np.zeros(shape=(self._n,), dtype=np.complex128),
)
)
else:
# non linear normalized case
self.K = K_non_lin_normalized
self.B = B_non_lin_normalized
self.C = C_non_lin_normalized
self.f_hierarchy = f_non_lin_hierarchy
self.f_hi_timed = f_non_lin_hierarchy_timed
self.x0 = np.concatenate(
(
self.get_x0_zeroTemp(),
np.zeros(shape=(self._n,), dtype=np.complex128),
)
)
self.f_const_args = F_Args_type(
eta=self.Eta,
O_vec=self.O_vec,
K=self.K,
B=self.B,
C=self.C,
M_down=self.M_down,
M_up=self.M_up,
spn=self.spn,
)
def __repr__(self):
return (
"sys_param : \n{}\n".format(self.SysP)
+ "int_param : \n{}\n".format(self.IntP)
+ "hi_param : \n{}\n".format(self.HiP)
+ "dimension of coupling Matrix : {}\n".format(self._dim)
+ "number of equations : {}\n".format(self.numEq)
)
def get_x0_zeroTemp(self):
if not len(self.SysP.psi0) == self._dim_sys:
raise RuntimeError(
"psiSys must match the dimension of the system Hilbertspace!"
)
res = np.zeros(self._dim, dtype=np.complex128)
res[: self._dim_sys] = self.SysP.psi0
return res
def get_args(self):
const_arg = {}
const_arg["z_num"] = 2 * self.Eta.get_num_y()
if self.EtaTherm is not None:
const_arg["stoc_temp_z_num"] = 2 * self.EtaTherm.get_num_y()
else:
const_arg["stoc_temp_z_num"] = 0
const_arg["t0"] = 0
const_arg["t1"] = self.IntP.t_max
const_arg["N"] = self.IntP.t_steps
const_arg["x0"] = self.x0
const_arg["f"] = self.f_hierarchy
const_arg["args"] = self.f_const_args
const_arg["integrator"] = self.IntP.integrator_name
if self.IntP.atol is not None:
const_arg["atol"] = self.IntP.atol
if self.IntP.rtol is not None:
const_arg["rtol"] = self.IntP.rtol
if self.IntP.order is not None:
const_arg["order"] = self.IntP.order
if self.IntP.nsteps is not None:
const_arg["nsteps"] = self.IntP.nsteps
if self.IntP.method is not None:
const_arg["method"] = self.IntP.method
if self.HiP.result_type == hid.RESULT_TYPE_ZEROTH_ORDER_ONLY:
const_arg["res_dim"] = (self._dim_sys,)
const_arg["x_to_res"] = partial(
x_to_res_zeroth_order_only, dim_sys=self._dim_sys
)
elif self.HiP.result_type == hid.RESULT_TYPE_ALL:
const_arg["res_dim"] = None
const_arg["x_to_res"] = None
elif self.HiP.result_type == hid.RESULT_TYPE_ZEROTH_ORDER_AND_ETA_LAMBDA:
if self.HiP.nonlinear is False:
raise RuntimeError("linear hierarchy has no ETA_LAMBDA")
const_arg["res_dim"] = (self._dim_sys + self._n,)
const_arg["x_to_res"] = partial(
x_to_res_zeroth_order_and_eta_lambda,
dim_sys=self._dim_sys,
num_bcf_terms=self._n,
)
else:
raise ValueError("unknown result_type '{}'".format(self.HiP.result_type))
const_arg["rand_skip"] = self.rand_skip
return const_arg
def measure_time_for_f_call(self, n=1):
const_arg = self.get_args()
arg = args_keyword_as_named_tupled(id=0, sub_seed=0)
args_dgl, kwargs = merge_arg_and_const_arg(arg, const_arg)
np.random.seed(arg.sub_seed)
rand_skip = kwargs["rand_skip"]
np.random.rand(rand_skip)
z_num = kwargs["z_num"]
stoc_temp_z_num = kwargs["stoc_temp_z_num"]
kwargs.pop("z_num")
kwargs.pop("stoc_temp_z_num")
kwargs.pop("sub_seed")
kwargs.pop("rand_skip")
z = uni_to_gauss(np.random.rand(z_num))
# this is a stoc proc instance, we set the process according to the RV z
args_dgl[0].new_process(z)
if stoc_temp_z_num != 0:
stoc_temp_z = uni_to_gauss(np.random.rand(stoc_temp_z_num))
args_dgl[-1].H_dyn[-1].new_process(stoc_temp_z)
del stoc_temp_z
f = self.f_hi_timed
x = self.x0
t0 = time()
for i in range(n):
t = np.random.rand() * self.IntP.t_max
dx_dt = f(t, x, *args_dgl)
print(dx_dt[:4])
t1 = time()
return (t1 - t0) / n
def integrate_simple(
self, data_name, data_path=".", overwrite=False, clear_pd=False
):
const_arg = self.get_args()
c_samples = progress.UnsignedIntValue()
c_int = progress.UnsignedIntValue()
m_samples = progress.UnsignedIntValue(val=self.number_of_samples)
m_int = progress.UnsignedIntValue(val=const_arg["N"])
md = hid.HIMetaData(hid_name=data_name, hid_path=data_path)
with md.get_HIData(key=self.hi_key) as data:
md.close()
if clear_pd:
log.info("clear HIData contained in {}/{}".format(data_path, data_name))
data.clear()
with progress.ProgressBarFancy(
count=[c_samples, c_int],
max_count=[m_samples, m_int],
prepend=["samples :", "integration :"],
interval=2,
) as p:
p.start()
skipped = 0
np.random.seed(self.HiP.seed)
np.random.rand(self.rand_skip)
subseeds_float = np.random.rand(self.number_of_samples) * 2 ** 32
for i in range(self.number_of_samples):
arg = args_keyword_as_named_tupled(
id=i, sub_seed=int(subseeds_float[i])
)
if overwrite:
log.info(
"calculation for idx {} TRIGGERED due to overwrite = True".format(
i
)
)
do_calculation = True
else:
do_calculation = not data.has_sample(idx=i)
if do_calculation:
log.info(
"idx {} NOT found in HiData, calculation TRIGGERED".format(
i
)
)
else:
log.info(
"idx {} found in HiData, calculation SKIPPED".format(i)
)
if do_calculation:
(t, psi_all, e_and_trb), z = _integrate_hierarchy(
arg, const_arg, c_int, m_int
)
if e_and_trb is not None:
print(e_and_trb[0])
print(e_and_trb[1])
raise NotImplementedError(
"no handling for incomplete integration data implemented"
)
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,
y=z,
)
else:
skipped += 1
c_samples.value = i + 1
p.reset(1)
log.info(
"{} out of {} arguments skipped!".format(
skipped, self.number_of_samples
)
)
def uni_to_gauss(x):
n = len(x) // 2
phi = x[:n] * 2 * np.pi
r = np.sqrt(-np.log(x[n:]))
return r * np.exp(1j * phi)
def complexRandGaussDistr(n, seed=None):
if seed != None:
np.random.seed(seed)
phi = np.random.rand(n) * 2 * np.pi
r = np.sqrt(-np.log(np.random.rand(n)))
z = r * np.exp(1j * phi)
# z = np.random.normal(size=2*n, scale=1/np.sqrt(2)).view(dtype=np.complex128)
return z
def complexSobolGaussDistr(n, seed=None):
if seed != None:
sbl.seed = seed
x = sbl.getNextSobol(dim=2 * n)
phi = x[::2] * 2 * np.pi
r = np.sqrt(-np.log(x[1::2]))
return r * np.exp(1j * phi)
Reference in a new issue View git blame Copy permalink