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) = / """ 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) = """ 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 = 1 with the properties = 0 = 0 = 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) stoc_temp_z = None 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) # t0, t1, N, f, args, x0, integrator, verbose, c, **kwargs return ode_wrapper.integrate_cplx(c=c, args=args_dgl, **kwargs), z, stoc_temp_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, temp_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, temp_y=temp_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)