diff --git a/pyproject.toml b/pyproject.toml index b888a3b..a462df2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,8 @@ dependencies = [ "scipy", "numba", "pytest", - "pytest-level" + "pytest-level", + "pytest-cov" ] requires-python = ">=3.12" authors = [ diff --git a/pytest.ini b/pytest.ini index be7a490..3f6d6c1 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,13 +1,15 @@ -[tool:pytest] +[pytest] testpaths = tests python_files = test_*.py python_classes = Test* python_functions = test_* +markers = + order: enforce execution order of tests addopts = --verbose + --level=1 --cov=mesohops --cov-report=term-missing --cov-report=html --cov-report=xml --cov-fail-under=70 - diff --git a/src/mesohops/basis/hops_aux.py b/src/mesohops/basis/hops_aux.py index 186b5a8..ab11e02 100644 --- a/src/mesohops/basis/hops_aux.py +++ b/src/mesohops/basis/hops_aux.py @@ -31,7 +31,6 @@ class AuxiliaryVector(Mapping): '_index', # Current index '__id_string', # Unique identifier string '__mode_digits', # Number of digits per mode - '__kmax', # Maximum hierarchy depth # --- Cached values --- '_sum', # Cached sum of vector elements diff --git a/src/mesohops/basis/hops_basis.py b/src/mesohops/basis/hops_basis.py index 3c094ce..2d6937c 100644 --- a/src/mesohops/basis/hops_basis.py +++ b/src/mesohops/basis/hops_basis.py @@ -935,7 +935,7 @@ def get_Z2_noise_sparse(self, z_step): if self.off_diagonal_couplings: # Get the noise associated with system-bath projection operators that # couple states in the current basis to a different state. - noise_t = (np.conj(z_step[0]) - 2j * np.real(z_step[1]))[ + noise_t = (np.conj(z_step[0]) - 1j * z_step[1])[ self.mode.list_rel_ind_off_diag_L2] # Get the noise memory drift associated with system-bath projection # operators that couple states in the current basis to a different state. diff --git a/src/mesohops/basis/hops_fluxfilters.py b/src/mesohops/basis/hops_fluxfilters.py index 4a322a0..110b5b2 100644 --- a/src/mesohops/basis/hops_fluxfilters.py +++ b/src/mesohops/basis/hops_fluxfilters.py @@ -332,6 +332,9 @@ def construct_filter_markov_up(self): filtered out) while False indicates otherwise (positioning is (mode, aux)). """ + if len(self.mode.list_absindex_mode) == 0: + return True + M2_mark_filtered_modes = np.array( [ np.array([param[m] for m in self.mode.list_absindex_mode]) @@ -389,6 +392,9 @@ def construct_filter_triangular_up(self): filtered out) while False indicates otherwise (positioning is (mode, aux)). """ + if len(self.mode.list_absindex_mode) == 0: + return True + M2_tri_filtered_modes = np.array( [ np.array([param[0][m] for m in self.mode.list_absindex_mode]) @@ -442,6 +448,9 @@ def construct_filter_longedge_up(self): filtered out) while False indicates otherwise (positioning is (mode, aux)). """ + if len(self.mode.list_absindex_mode) == 0: + return True + M2_le_filtered_modes = np.array( [ np.array([param[0][m] for m in self.mode.list_absindex_mode]) diff --git a/src/mesohops/basis/hops_hierarchy.py b/src/mesohops/basis/hops_hierarchy.py index 9b92c99..0709d13 100644 --- a/src/mesohops/basis/hops_hierarchy.py +++ b/src/mesohops/basis/hops_hierarchy.py @@ -3,7 +3,8 @@ import numpy as np -from mesohops.basis.hierarchy_functions import (filter_aux_longedge, filter_aux_triangular, filter_markovian) +from mesohops.basis.hierarchy_functions import (filter_aux_longedge, + filter_aux_triangular, filter_markovian) from mesohops.basis.hops_aux import AuxiliaryVector as AuxVec from mesohops.util.dynamic_dict import Dict_wDefaults from mesohops.util.exceptions import AuxError, UnsupportedRequest @@ -12,7 +13,6 @@ __author__ = "D. I. G. Bennett, L. Varvelo, J. K. Lynd" __version__ = "1.2" - HIERARCHY_DICT_DEFAULT = {"MAXHIER": int(3), "TERMINATOR": False, "STATIC_FILTERS": []} HIERARCHY_DICT_TYPES = dict( @@ -31,33 +31,33 @@ class HopsHierarchy(Dict_wDefaults): __slots__ = ( # --- Core basis components --- - 'system', # System parameters and operators (HopsSystem) - 'n_hmodes', # Number of hierarchy modes - '__ndim', # System dimension (for internal use) + 'system', # System parameters and operators (HopsSystem) + 'n_hmodes', # Number of hierarchy modes # --- Auxiliary list management --- - '_auxiliary_list', # List of current auxiliary vectors (main storage) - '__previous_auxiliary_list', # Auxiliary list from previous step - '__list_aux_stable', # Auxiliaries stable between steps - '__list_aux_remove', # Auxiliaries to remove in update - '__list_aux_add', # Auxiliaries to add in update - '__previous_list_auxstable_index', # Indices of stable auxiliaries from previous step - '_previous_list_modes_in_use', # Modes in use from previous step + '_auxiliary_list', # List of current auxiliary vectors (main storage) + '__previous_auxiliary_list', # Auxiliary list from previous step + '__list_aux_stable', # Auxiliaries stable between steps + '__list_aux_remove', # Auxiliaries to remove in update + '__list_aux_add', # Auxiliaries to add in update + '__previous_list_auxstable_index', + # Indices of stable auxiliaries from previous step + '_previous_list_modes_in_use', # Modes in use from previous step # --- Filter management flags --- - 'only_markovian_filter', # True if only Markovian filter is used, else False + 'only_markovian_filter', # True if only Markovian filter is used, else False # --- Parameter management --- - '_default_param', # Default parameter dictionary - '_param_types', # Parameter type dictionary + '_default_param', # Default parameter dictionary + '_param_types', # Parameter type dictionary # --- Connection and indexing dictionaries --- - '_new_aux_index_conn_by_mode', # New auxiliary index connections by mode - '_new_aux_id_conn_by_mode', # New auxiliary ID connections by mode - '_stable_aux_id_conn_by_mode', # Stable auxiliary ID connections by mode - '_dict_aux_by_id', # Dictionary mapping auxiliary IDs to objects - '_list_modes_in_use', # List of modes currently in use - '_count_by_modes', # Count of auxiliaries by mode + '_new_aux_index_conn_by_mode', # New auxiliary index connections by mode + '_new_aux_id_conn_by_mode', # New auxiliary ID connections by mode + '_stable_aux_id_conn_by_mode', # Stable auxiliary ID connections by mode + '_dict_aux_by_id', # Dictionary mapping auxiliary IDs to objects + '_list_modes_in_use', # List of modes currently in use + '_count_by_modes', # Count of auxiliaries by mode ) def __init__(self, hierarchy_param, system_param): @@ -106,13 +106,43 @@ def __init__(self, hierarchy_param, system_param): self.param["MAXHIER"] = 255 self.n_hmodes = system_param["N_HMODES"] self._auxiliary_list = [] - self._new_aux_index_conn_by_mode = {mode: {} for mode in range(system_param["N_HMODES"])} - self._new_aux_id_conn_by_mode = {mode: {} for mode in range(system_param["N_HMODES"])} - self._stable_aux_id_conn_by_mode = {mode: {} for mode in range(system_param["N_HMODES"])} + self._new_aux_index_conn_by_mode = {mode: {} for mode in + range(system_param["N_HMODES"])} + self._new_aux_id_conn_by_mode = {mode: {} for mode in + range(system_param["N_HMODES"])} + self._stable_aux_id_conn_by_mode = {mode: {} for mode in + range(system_param["N_HMODES"])} self._dict_aux_by_id = {} self._list_modes_in_use = [] self._count_by_modes = {} + # Prevent the application of static filters with the wrong number of modes. + if "STATIC_FILTERS" in self.param.keys(): + for (filter_i, param_i) in self.param["STATIC_FILTERS"]: + if filter_i == "Triangular": + if not len(param_i[0]) == self.n_hmodes: + raise UnsupportedRequest("The number of entries in the list of " + "static filter booleans does not match the " + "number of hierarchy modes", + "hierarchy_static_filter", + override=True) + elif filter_i == "LongEdge": + if not len(param_i[0]) == self.n_hmodes: + raise UnsupportedRequest("The number of entries in the list of " + "static filter booleans does not match the " + "number of hierarchy modes", + "hierarchy_static_filter", + override=True) + elif filter_i == "Markovian": + if not len(param_i) == self.n_hmodes: + raise UnsupportedRequest("The number of entries in the list of " + "static filter booleans does not match the " + "number of hierarchy modes", + "hierarchy_static_filter", + override=True) + else: + raise UnsupportedRequest(filter_i, "hierarchy_static_filter") + def initialize(self, flag_adaptive): """ Initializes the hierarchy. @@ -121,7 +151,7 @@ def initialize(self, flag_adaptive): ---------- 1. flag_adaptive : bool True indicates an adaptive calculation while False indicates - otherwise. + otherwise. Returns ------- @@ -165,7 +195,8 @@ def initialize(self, flag_adaptive): # Nsteps). At the moment, though, I'm skipping this and allowing # the hierarchy to control its own growth. self.auxiliary_list = [AuxVec([], self.n_hmodes)] - if np.any([name != "Markovian" for (name, param) in self.param["STATIC_FILTERS"]]): + if np.any([name != "Markovian" for (name, param) in + self.param["STATIC_FILTERS"]]): self.only_markovian_filter = False else: self.only_markovian_filter = True @@ -248,8 +279,8 @@ def apply_filter(self, list_aux, filter_name, params): # Update STATIC_FILTERS parameters if needed # ------------------------------------------ if not ( - [filter_name, params] in self.param["STATIC_FILTERS"] - or (filter_name, params) in self.param["STATIC_FILTERS"] + [filter_name, params] in self.param["STATIC_FILTERS"] + or (filter_name, params) in self.param["STATIC_FILTERS"] ): self.param["STATIC_FILTERS"].append((filter_name, params)) @@ -260,11 +291,11 @@ def _aux_index(self, aux): Returns the index value for a given auxiliary. The important thing is that this function is aware of whether the calculation should be using an 'absolute' index or a 'relative' index. - + Absolute index: no matter which auxiliaries are in the hierarchy, the index value does not change. This is a useful approach when trying - to do things like dynamic filtering. - + to do things like dynamic filtering. + Relative index: This is the more intuitive indexing scheme which only keeps track of the auxiliary vectors that are actually in the hierarchy. @@ -287,7 +318,7 @@ def _aux_index(self, aux): def _const_aux_edge(absindex_mode, depth, n_hmodes): """ Creates an auxiliary object for an edge node at - a particular depth along a given mode. + a particular depth along a given mode. Parameters ---------- @@ -365,7 +396,7 @@ def define_markovian_filtered_triangular_hierarchy(n_hmodes, maxhier, list_aux = [] # Loop over hierarchy depths at which the Markovian filter does not apply - for k in [0,1]: + for k in [0, 1]: for aux_raw in it.combinations_with_replacement(np.arange(n_hmodes), k): count = Counter(aux_raw) list_aux.append( @@ -375,7 +406,7 @@ def define_markovian_filtered_triangular_hierarchy(n_hmodes, maxhier, # Generate an array of exclusively the non-Markovian-filtered modes M1_modes_filtered = np.arange(n_hmodes)[list_not_boolean_mark] # Loop over hierarchy depths at which the Markovian filter applies - for k in np.arange(2, maxhier+1): + for k in np.arange(2, maxhier + 1): # At each depth, add to list_aux all possible combinations of only the # non-Markovian-filtered modes for aux_raw in it.combinations_with_replacement(M1_modes_filtered, k): @@ -402,7 +433,7 @@ def __update_count(self, aux, type): None """ if type == 'add': - self._count_by_modes |= {mode:(self._count_by_modes[mode] + 1) + self._count_by_modes |= {mode: (self._count_by_modes[mode] + 1) if mode in self._count_by_modes.keys() else 1 for mode in aux.keys()} elif type == 'remove': self._count_by_modes |= {mode: (self._count_by_modes[mode] - 1) @@ -428,7 +459,7 @@ def __update_modes_in_use(self): for (mode, value) in self._count_by_modes.items(): if value == 0: list_to_remove.append(mode) - elif value>0: + elif value > 0: list_modes_in_use.append(mode) else: print(f'ERROR: _count_by_modes is negative for mode {mode}') @@ -449,24 +480,30 @@ def add_connections(self): None """ for mode in self._previous_list_modes_in_use: - self._stable_aux_id_conn_by_mode[mode].update(self._new_aux_id_conn_by_mode[mode]) - self._new_aux_index_conn_by_mode = {mode: {} for mode in self._list_modes_in_use} + self._stable_aux_id_conn_by_mode[mode].update( + self._new_aux_id_conn_by_mode[mode]) + self._new_aux_index_conn_by_mode = {mode: {} for mode in + self._list_modes_in_use} self._new_aux_id_conn_by_mode = {mode: {} for mode in self._list_modes_in_use} for aux in self.list_aux_add: sum_aux = np.sum(aux) # add connections to k+1 if sum_aux < self.param['MAXHIER']: - list_id_p1, list_value_connects_p1, list_mode_connects_p1 = aux.get_list_id_up(self._list_modes_in_use) - for (rel_ind,my_id) in enumerate(list_id_p1): + list_id_p1, list_value_connects_p1, list_mode_connects_p1 = aux.get_list_id_up( + self._list_modes_in_use) + for (rel_ind, my_id) in enumerate(list_id_p1): try: aux_p1 = self._dict_aux_by_id[my_id] - aux.add_aux_connect(list_mode_connects_p1[rel_ind],aux_p1,1) - - #We simply keep track of the index connections of new auxiliaries - #Note: It does not matter that the indices will change because we only use this dictionary - #once, immediately after it is created in eom.ksuper - self._new_aux_id_conn_by_mode[list_mode_connects_p1[rel_ind]][aux.id] = [aux_p1.id,list_value_connects_p1[rel_ind] + 1] - self._new_aux_index_conn_by_mode[list_mode_connects_p1[rel_ind]][aux._index] = [aux_p1._index,list_value_connects_p1[rel_ind] + 1] + aux.add_aux_connect(list_mode_connects_p1[rel_ind], aux_p1, 1) + + # We simply keep track of the index connections of new auxiliaries + # Note: It does not matter that the indices will change because we only use this dictionary + # once, immediately after it is created in eom.ksuper + self._new_aux_id_conn_by_mode[list_mode_connects_p1[rel_ind]][ + aux.id] = [aux_p1.id, list_value_connects_p1[rel_ind] + 1] + self._new_aux_index_conn_by_mode[ + list_mode_connects_p1[rel_ind]][aux._index] = [ + aux_p1._index, list_value_connects_p1[rel_ind] + 1] except: pass @@ -474,23 +511,25 @@ def add_connections(self): if sum_aux > 0: list_id_m1, list_value_connects_m1, list_mode_connects_m1 = aux.get_list_id_down() - for (index,id_m1) in enumerate(list_id_m1): + for (index, id_m1) in enumerate(list_id_m1): try: aux_m1 = self._dict_aux_by_id[id_m1] aux.add_aux_connect(list_mode_connects_m1[index], aux_m1, -1) - self._new_aux_id_conn_by_mode[list_mode_connects_m1[index]][aux_m1.id] = [aux.id,list_value_connects_m1[index]] - self._new_aux_index_conn_by_mode[list_mode_connects_m1[index]][aux_m1._index] = [aux._index,list_value_connects_m1[index]] + self._new_aux_id_conn_by_mode[list_mode_connects_m1[index]][ + aux_m1.id] = [aux.id, list_value_connects_m1[index]] + self._new_aux_index_conn_by_mode[list_mode_connects_m1[index]][ + aux_m1._index] = [aux._index, list_value_connects_m1[index]] except: pass - + @property def new_aux_index_conn_by_mode(self): return self._new_aux_index_conn_by_mode - + @property def stable_aux_id_conn_by_mode(self): - return self._stable_aux_id_conn_by_mode - + return self._stable_aux_id_conn_by_mode + @property def size(self): return len(self.auxiliary_list) @@ -520,14 +559,15 @@ def auxiliary_list(self, aux_list): self.__list_aux_remove = list_aux_remove self.__list_aux_add = list(set_aux_add) self.__list_aux_add.sort() - + set_aux_stable = set(self.auxiliary_list).intersection(set(aux_list)) - self.__list_aux_stable = [aux for aux in self.auxiliary_list if aux in set_aux_stable] - self.__previous_list_auxstable_index = [aux._index for aux in self.__list_aux_stable] - + self.__list_aux_stable = [aux for aux in self.auxiliary_list if + aux in set_aux_stable] + self.__previous_list_auxstable_index = [aux._index for aux in + self.__list_aux_stable] + self._previous_list_modes_in_use = self._list_modes_in_use.copy() - - + if set(aux_list) != set(self.__previous_auxiliary_list): # Prepare New Auxiliary List # -------------------------- @@ -543,7 +583,7 @@ def auxiliary_list(self, aux_list): self._auxiliary_list = aux_list if not aux_list[0].id == '': raise AuxError("Zero Vector not contained in list_aux!") - + # Update modes_in_use self.__update_modes_in_use() @@ -551,29 +591,27 @@ def auxiliary_list(self, aux_list): for (index, aux) in enumerate(aux_list): aux._index = index - - # Add auxiliary connections self.add_connections() - + # Delete stable connections of deleted auxiliaries - + for aux in list_aux_remove: - for (mode,aux_p1) in aux.dict_aux_p1.items(): + for (mode, aux_p1) in aux.dict_aux_p1.items(): try: self._stable_aux_id_conn_by_mode[mode].pop(aux.id) except: pass - for (mode,aux_m1) in aux.dict_aux_m1.items(): + for (mode, aux_m1) in aux.dict_aux_m1.items(): try: self._stable_aux_id_conn_by_mode[mode].pop(aux_m1.id) except: pass - + # Remove auxiliary connections for aux in set_aux_remove: aux.remove_pointers() - + @property def list_aux_stable(self): return self.__list_aux_stable diff --git a/src/mesohops/basis/hops_system.py b/src/mesohops/basis/hops_system.py index 83b58d5..50f5d81 100644 --- a/src/mesohops/basis/hops_system.py +++ b/src/mesohops/basis/hops_system.py @@ -1,10 +1,15 @@ -from collections import Counter +from __future__ import annotations + +import os +import pickle +from pathlib import Path +from typing import Any, Sequence import numpy as np import scipy as sp from scipy import sparse -from mesohops.util.helper_functions import array_to_tuple +from mesohops.basis.system_functions import initialize_system_dict __title__ = "System Class" __author__ = "D. I. G. Bennett, L. Varvelo, J. K. Lynd, B. Z. Citty" @@ -29,7 +34,6 @@ class HopsSystem: 'adaptive', # Adaptive flag (True if adaptive basis is used) '__list_add_state', # States to add in update '__list_stable_state', # States stable between updates - '_list_off_diag', # Off-diagonal elements of the Hamiltonian '_list_boundary_state', # States coupled to basis by Hamiltonian # --- Indexing of modes & L-operators in the current basis --- @@ -40,13 +44,14 @@ class HopsSystem: '__dict_relindex_states', # Relative state indices ) - def __init__(self, system_param): + def __init__(self, system_param: dict[str, Any] | str | os.PathLike[str] | Path) -> None: """ Inputs ------ - 1. system_param : dict - Dictionary with the system and system-bath coupling - parameters defined. + 1. system_param : dict | str | Path + Either a dictionary with the system and system-bath coupling + parameters defined or a str/pathlike that points to a file + that has been generated with the method save_dict_param(). a. HAMILTONIAN : np.array Array that contains the system Hamiltonian. b. GW_SYSBATH : list(complex) @@ -120,176 +125,17 @@ def __init__(self, system_param): that will get its own super-operators, but since we have no use-case yet this has not been implemented. """ - self.param = self._initialize_system_dict(system_param) + if isinstance(system_param, (str, os.PathLike)): + self.param = pickle.load(open(system_param, "rb")) + elif isinstance(system_param, dict): + self.param = initialize_system_dict(system_param) + else: + raise TypeError("system_param must be a dictionary or a file path.") self.__ndim = self.param["NSTATES"] self.__previous_state_list = None self.__state_list = [] - def _initialize_system_dict(self, system_param): - """ - Extends the user input to the complete set of parameters defined above. - - Parameters - ---------- - 1. system_param : dict - Dictionary with the system and system-bath coupling - parameters defined. - - Returns - ------- - 1. param_dict : dict - Dictionary containing the user input and the derived parameters. - """ - param_dict = system_param - if(sparse.issparse(system_param["HAMILTONIAN"])): - param_dict["NSTATES"] = sparse.coo_matrix.get_shape(system_param["HAMILTONIAN"])[0] - else: - param_dict["NSTATES"] = len(system_param["HAMILTONIAN"][0]) - param_dict["N_HMODES"] = len(system_param["GW_SYSBATH"]) - param_dict["G"] = np.array([g for (g, w) in system_param["GW_SYSBATH"]]) - param_dict["W"] = np.array([w for (g, w) in system_param["GW_SYSBATH"]]) - param_dict["LIST_STATE_INDICES_BY_HMODE"] = [ - self._get_states_from_L2(L2) for L2 in param_dict["L_HIER"] - ] - param_dict["LIST_HMODE_INDICES_BY_STATE"] = [[] for i in range(param_dict["NSTATES"])] - for (hmode,state_indices) in enumerate(param_dict["LIST_STATE_INDICES_BY_HMODE"]): - for state in state_indices: - param_dict["LIST_HMODE_INDICES_BY_STATE"][state].append(hmode) - - - param_dict["SPARSE_HAMILTONIAN"] = sparse.csc_array(param_dict["HAMILTONIAN"]) - param_dict["SPARSE_HAMILTONIAN"].eliminate_zeros() - - sparse_ham = sparse.coo_matrix(system_param["HAMILTONIAN"]) - param_dict["COUPLED_STATES"] = [[] for state in range(param_dict["NSTATES"])] - for i in range(len(sparse_ham.row)): - param_dict["COUPLED_STATES"][sparse_ham.row[i]].append(sparse_ham.col[i]) - for i in range(param_dict["NSTATES"]): - param_dict["COUPLED_STATES"][i] = sorted(param_dict["COUPLED_STATES"][i]) - - # Checks for low-temperature correction terms - if there are none, initialize - # empty lists as placeholders: - if not "L_LT_CORR" in param_dict.keys(): - param_dict["L_LT_CORR"] = [] - param_dict["PARAM_LT_CORR"] = [] - - # Define the Hierarchy Operator Values - # ------------------------------------ - # Since arrays and lists are not hashable, we will turn our operators - # into tuples in order to conveniently define a number of indexing - # parameters. - - # Creates list of unique l2 tuples in order they appear in "L_HIER" - l2_as_tuples = [array_to_tuple(L2) for L2 in param_dict["L_HIER"]] - list_unique_l2_as_tuples = list(Counter(l2_as_tuples)) - param_dict["N_L2"] = len(set(list_unique_l2_as_tuples)) - # Generates a list of destination state indices linked to each state index by - # one of the L-operators. - param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"] = [[] for s in range( - param_dict["NSTATES"])] - for l in list_unique_l2_as_tuples: - for j in range(len(l[0])): - param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"][l[0][j]].append( - l[1][j]) - param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"] = [list(set( - list_dest_state)) for list_dest_state in - param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"]] - - # Creates L2 indexing parameters - param_dict["LIST_INDEX_L2_BY_HMODE"] = [ - None for i in range(param_dict["N_HMODES"]) - ] - param_dict["LIST_L2_COO"] = [0]*param_dict["N_L2"] - param_dict["LIST_LT_PARAM"] = [0]*param_dict["N_L2"] - param_dict["LIST_HMODE_BY_INDEX_L2"] = [[]]*param_dict["N_L2"] - - param_dict["LIST_STATE_INDICES_BY_INDEX_L2"] = [] - dict_unique_l2 = {} - unique_l2_insertion = 0 - for (i,l_2) in enumerate(l2_as_tuples): - try: - index_l2_unique = dict_unique_l2[l_2] - param_dict["LIST_INDEX_L2_BY_HMODE"][i] = index_l2_unique - param_dict["LIST_HMODE_BY_INDEX_L2"][index_l2_unique].append(i) - - except: - dict_unique_l2[l_2] = unique_l2_insertion - param_dict["LIST_INDEX_L2_BY_HMODE"][i] = unique_l2_insertion - param_dict["LIST_HMODE_BY_INDEX_L2"][unique_l2_insertion].append(i) - tmp = sp.sparse.coo_matrix(param_dict["L_HIER"][i]) - tmp.eliminate_zeros() - param_dict["LIST_L2_COO"][unique_l2_insertion] = tmp - param_dict["LIST_STATE_INDICES_BY_INDEX_L2"].append( - param_dict["LIST_STATE_INDICES_BY_HMODE"][i] - ) - unique_l2_insertion += 1 - - self._list_off_diag = np.array([not np.allclose(L2.col, L2.row) - for L2 in param_dict["LIST_L2_COO"]]) - - param_dict["LIST_INDEX_L2_BY_STATE_INDICES"] = [[] for i in range(param_dict["NSTATES"])] - for (index_L2,state_indices) in enumerate(param_dict["LIST_STATE_INDICES_BY_INDEX_L2"]): - for state in state_indices: - param_dict["LIST_INDEX_L2_BY_STATE_INDICES"][state].append(index_L2) - - l2_LT_CORR_as_tuples = [array_to_tuple(l) for l in - param_dict["L_LT_CORR"]] - dict_index_unique_l2_as_tuples = {} - for (i,l) in enumerate(list_unique_l2_as_tuples): - dict_index_unique_l2_as_tuples[l] = i - - # Build a list of low-temperature coefficients guaranteed to be in the same - # order as the associated unique sparse L2 operators. - for (i,l) in enumerate(l2_LT_CORR_as_tuples): - try: - mister_mc_index = dict_index_unique_l2_as_tuples[l] - param_dict["LIST_LT_PARAM"][mister_mc_index] += param_dict["PARAM_LT_CORR"][i] - except: - print("WARNING: the list of low-temperature correction " - "L-operators contains an L-operator not associated with any " - "existing thermal environment. This low-temperature " - "correction factor will be discarded!") - - # Define the Noise1 Operator Values - # --------------------------------- - param_dict["LIST_NMODE1_BY_INDEX_L2"] = [[] for i in range(len(l2_as_tuples))] - param_dict["LIST_INDEX_L2_BY_NMODE1"] = [ - None for i in range(len(param_dict["PARAM_NOISE1"])) - ] - l2_as_tuples = [array_to_tuple(l) for l in param_dict["L_NOISE1"]] - for (i, l_2) in enumerate(l2_as_tuples): - try: - index_l2_unique = dict_unique_l2[l_2] - param_dict["LIST_INDEX_L2_BY_NMODE1"][i] = index_l2_unique - param_dict["LIST_NMODE1_BY_INDEX_L2"][index_l2_unique].append(i) - except: - print("WARNING: the list of noise 1 L-operators contains an L-operator " - "not associated with any existing thermal environment. The noise " - "associated with this L-operator will be discarded!") - - # Define the Noise2 Operator Values - # --------------------------------- - if "L_NOISE2" in param_dict.keys(): - param_dict["LIST_NMODE2_BY_INDEX_L2"] = [[] for i in range(len(l2_as_tuples))] - param_dict["LIST_INDEX_L2_BY_NMODE2"] = [ - None for i in range(len(param_dict["PARAM_NOISE2"])) - ] - l2_as_tuples = [array_to_tuple(l) for l in param_dict["L_NOISE2"]] - - for (i, l_2) in enumerate(l2_as_tuples): - - try: - index_l2_unique = dict_unique_l2[l_2] - param_dict["LIST_INDEX_L2_BY_NMODE2"][i] = index_l2_unique - param_dict["LIST_NMODE2_BY_INDEX_L2"][index_l2_unique].append(i) - except: - print("WARNING: the list of noise 2 L-operators contains an L-operator " - "not associated with any existing thermal environment. The noise " - "associated with this L-operator will be discarded!") - - return param_dict - - def initialize(self, flag_adaptive, psi_0): + def initialize(self, flag_adaptive: bool, psi_0: np.ndarray) -> None: """ Creates a state list depending on whether the calculation is adaptive or not. @@ -313,51 +159,52 @@ def initialize(self, flag_adaptive, psi_0): else: self.state_list = np.arange(self.__ndim) - @staticmethod - def _get_states_from_L2(lop): + def save_dict_param(self, filepath: str | os.PathLike[str] | Path) -> None: """ - Fetches the states that the L operators interacts with. + Serialize the system parameters to a file. + + This method saves the current system parameters stored in `self.param` to the specified + file path using the pickle serialization format. This allows for easy storage and retrieval + of the system's configuration. Parameters ---------- - 1. lop : np.array(complex) - L2 operator. + 1. filepath : str or os.PathLike + The path to the file where the system parameters will be saved. Returns ------- - 1. tuple : tuple - Tuple of states that correspond to the specific L operator. + None """ - - i_x, i_y = np.nonzero(lop) - return tuple(set(i_x) | set(i_y)) + with open(filepath, "wb") as f: + pickle.dump(self.param, f) @property - def size(self): + def size(self) -> int: return len(self.__state_list) @property - def state_list(self): + def state_list(self) -> np.ndarray | list: return self.__state_list @property - def list_destination_state(self): + def list_destination_state(self) -> np.ndarray: return self.__list_destination_state @property - def list_boundary_state(self): + def list_boundary_state(self) -> list[int]: return self._list_boundary_state @property - def list_sc(self): + def list_sc(self) -> list[int]: list_boundary_lop = list(set(self.list_destination_state) - set(self.state_list)) return list(set(list_boundary_lop) | set(self.list_boundary_state)) @property - def dict_relative_index_by_state(self): + def dict_relative_index_by_state(self) -> dict[int, int]: return self.__dict_relindex_states @state_list.setter - def state_list(self, new_state_list): + def state_list(self, new_state_list: Sequence[int] | np.ndarray) -> None: # Construct information about previous timestep # -------------------------------------------- self.__previous_state_list = self.__state_list @@ -439,40 +286,42 @@ def state_list(self, new_state_list): self._list_boundary_state = list(set([state_conn for conn_list in self._list_boundary_state for state_conn in conn_list ]) - set(self.state_list)) @property - def previous_state_list(self): + def previous_state_list(self) -> np.ndarray | None: return self.__previous_state_list @property - def list_stable_state(self): + def list_stable_state(self) -> np.ndarray | list: return self.__list_stable_state @property - def list_add_state(self): + def list_add_state(self) -> np.ndarray | list: return self.__list_add_state @property - def hamiltonian(self): + def hamiltonian(self) -> sp.sparse.spmatrix | np.ndarray: return self._hamiltonian @property - def list_absindex_state_modes(self): + def list_absindex_state_modes(self) -> np.ndarray: return self.__list_absindex_state_modes @property - def list_absindex_new_state_modes(self): + def list_absindex_new_state_modes(self) -> np.ndarray: return self.__list_absindex_new_state_modes @property - def list_absindex_L2_active(self): + def list_absindex_L2_active(self) -> np.ndarray: return self.__list_absindex_L2_active @property - def list_lt_corr_param(self): + def list_lt_corr_param(self) -> np.ndarray: return self._list_lt_corr_param @property - def list_off_diag(self): - return self._list_off_diag + def list_off_diag(self) -> np.ndarray: + return self.param["list_L2_off_diag"] @staticmethod - def reduce_sparse_matrix(coo_mat, state_list): + def reduce_sparse_matrix( + coo_mat: sp.sparse.spmatrix, state_list: Sequence[int] + ) -> sp.sparse.coo_matrix: """ Takes in a sparse matrix and list which represents the absolute state to a new relative state represented in a sparse matrix. diff --git a/src/mesohops/basis/system_functions.py b/src/mesohops/basis/system_functions.py new file mode 100644 index 0000000..f9f4bbb --- /dev/null +++ b/src/mesohops/basis/system_functions.py @@ -0,0 +1,172 @@ +from __future__ import annotations + +from collections import Counter +from typing import Any, Dict + +import numpy as np +from scipy import sparse +from mesohops.util.helper_functions import get_states_from_L2, array_to_tuple + +def initialize_system_dict(system_param: Dict[str, Any]) -> Dict[str, Any]: + """ + Extends the user input to the complete set of parameters defined above. + + Parameters + ---------- + 1. system_param : dict + Dictionary with the system and system-bath coupling + parameters defined. + + Returns + ------- + 1. param_dict : dict + Dictionary containing the user input and the derived parameters. + """ + param_dict = system_param + if(sparse.issparse(system_param["HAMILTONIAN"])): + param_dict["NSTATES"] = sparse.coo_matrix.get_shape(system_param["HAMILTONIAN"])[0] + else: + param_dict["NSTATES"] = len(system_param["HAMILTONIAN"][0]) + param_dict["N_HMODES"] = len(system_param["GW_SYSBATH"]) + param_dict["G"] = np.array([g for (g, w) in system_param["GW_SYSBATH"]]) + param_dict["W"] = np.array([w for (g, w) in system_param["GW_SYSBATH"]]) + param_dict["LIST_STATE_INDICES_BY_HMODE"] = [ + get_states_from_L2(L2) for L2 in param_dict["L_HIER"] + ] + param_dict["LIST_HMODE_INDICES_BY_STATE"] = [[] for i in range(param_dict["NSTATES"])] + for (hmode ,state_indices) in enumerate(param_dict["LIST_STATE_INDICES_BY_HMODE"]): + for state in state_indices: + param_dict["LIST_HMODE_INDICES_BY_STATE"][state].append(hmode) + + + param_dict["SPARSE_HAMILTONIAN"] = sparse.csc_array(param_dict["HAMILTONIAN"]) + param_dict["SPARSE_HAMILTONIAN"].eliminate_zeros() + + sparse_ham = sparse.coo_matrix(system_param["HAMILTONIAN"]) + param_dict["COUPLED_STATES"] = [[] for state in range(param_dict["NSTATES"])] + for i in range(len(sparse_ham.row)): + param_dict["COUPLED_STATES"][sparse_ham.row[i]].append(sparse_ham.col[i]) + for i in range(param_dict["NSTATES"]): + param_dict["COUPLED_STATES"][i] = sorted(param_dict["COUPLED_STATES"][i]) + + # Checks for low-temperature correction terms - if there are none, initialize + # empty lists as placeholders: + if not "L_LT_CORR" in param_dict.keys(): + param_dict["L_LT_CORR"] = [] + param_dict["PARAM_LT_CORR"] = [] + + # Define the Hierarchy Operator Values + # ------------------------------------ + # Since arrays and lists are not hashable, we will turn our operators + # into tuples in order to conveniently define a number of indexing + # parameters. + + # Creates list of unique l2 tuples in order they appear in "L_HIER" + l2_as_tuples = [array_to_tuple(L2) for L2 in param_dict["L_HIER"]] + list_unique_l2_as_tuples = list(Counter(l2_as_tuples)) + param_dict["N_L2"] = len(set(list_unique_l2_as_tuples)) + # Generates a list of destination state indices linked to each state index by + # one of the L-operators. + param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"] = [[] for s in range( + param_dict["NSTATES"])] + for l in list_unique_l2_as_tuples: + for j in range(len(l[0])): + param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"][l[0][j]].append( + l[1][j]) + param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"] = [list(set( + list_dest_state)) for list_dest_state in + param_dict["LIST_DESTINATION_STATES_BY_STATE_INDEX"]] + + # Creates L2 indexing parameters + param_dict["LIST_INDEX_L2_BY_HMODE"] = [ + None for i in range(param_dict["N_HMODES"]) + ] + param_dict["LIST_L2_COO"] = [0 ] *param_dict["N_L2"] + param_dict["LIST_LT_PARAM"] = [0 ] *param_dict["N_L2"] + param_dict["LIST_HMODE_BY_INDEX_L2"] = [[] ] *param_dict["N_L2"] + + param_dict["LIST_STATE_INDICES_BY_INDEX_L2"] = [] + dict_unique_l2 = {} + unique_l2_insertion = 0 + for (i ,l_2) in enumerate(l2_as_tuples): + try: + index_l2_unique = dict_unique_l2[l_2] + param_dict["LIST_INDEX_L2_BY_HMODE"][i] = index_l2_unique + param_dict["LIST_HMODE_BY_INDEX_L2"][index_l2_unique].append(i) + + except: + dict_unique_l2[l_2] = unique_l2_insertion + param_dict["LIST_INDEX_L2_BY_HMODE"][i] = unique_l2_insertion + param_dict["LIST_HMODE_BY_INDEX_L2"][unique_l2_insertion].append(i) + tmp = sparse.coo_matrix(param_dict["L_HIER"][i]) + tmp.eliminate_zeros() + param_dict["LIST_L2_COO"][unique_l2_insertion] = tmp + param_dict["LIST_STATE_INDICES_BY_INDEX_L2"].append( + param_dict["LIST_STATE_INDICES_BY_HMODE"][i] + ) + unique_l2_insertion += 1 + + param_dict["list_L2_off_diag"] = np.array([not np.allclose(L2.col, L2.row) + for L2 in param_dict["LIST_L2_COO"]]) + + param_dict["LIST_INDEX_L2_BY_STATE_INDICES"] = [[] for i in range(param_dict["NSTATES"])] + for (index_L2 ,state_indices) in enumerate(param_dict["LIST_STATE_INDICES_BY_INDEX_L2"]): + for state in state_indices: + param_dict["LIST_INDEX_L2_BY_STATE_INDICES"][state].append(index_L2) + + l2_LT_CORR_as_tuples = [array_to_tuple(l) for l in + param_dict["L_LT_CORR"]] + dict_index_unique_l2_as_tuples = {} + for (i ,l) in enumerate(list_unique_l2_as_tuples): + dict_index_unique_l2_as_tuples[l] = i + + # Build a list of low-temperature coefficients guaranteed to be in the same + # order as the associated unique sparse L2 operators. + for (i ,l) in enumerate(l2_LT_CORR_as_tuples): + try: + mister_mc_index = dict_index_unique_l2_as_tuples[l] + param_dict["LIST_LT_PARAM"][mister_mc_index] += param_dict["PARAM_LT_CORR"][i] + except: + print("WARNING: the list of low-temperature correction " + "L-operators contains an L-operator not associated with any " + "existing thermal environment. This low-temperature " + "correction factor will be discarded!") + + # Define the Noise1 Operator Values + # --------------------------------- + param_dict["LIST_NMODE1_BY_INDEX_L2"] = [[] for i in range(len(l2_as_tuples))] + param_dict["LIST_INDEX_L2_BY_NMODE1"] = [ + None for i in range(len(param_dict["PARAM_NOISE1"])) + ] + l2_as_tuples = [array_to_tuple(l) for l in param_dict["L_NOISE1"]] + for (i, l_2) in enumerate(l2_as_tuples): + try: + index_l2_unique = dict_unique_l2[l_2] + param_dict["LIST_INDEX_L2_BY_NMODE1"][i] = index_l2_unique + param_dict["LIST_NMODE1_BY_INDEX_L2"][index_l2_unique].append(i) + except: + print("WARNING: the list of noise 1 L-operators contains an L-operator " + "not associated with any existing thermal environment. The noise " + "associated with this L-operator will be discarded!") + + # Define the Noise2 Operator Values + # --------------------------------- + if "L_NOISE2" in param_dict.keys(): + param_dict["LIST_NMODE2_BY_INDEX_L2"] = [[] for i in range(len(l2_as_tuples))] + param_dict["LIST_INDEX_L2_BY_NMODE2"] = [ + None for i in range(len(param_dict["PARAM_NOISE2"])) + ] + l2_as_tuples = [array_to_tuple(l) for l in param_dict["L_NOISE2"]] + + for (i, l_2) in enumerate(l2_as_tuples): + + try: + index_l2_unique = dict_unique_l2[l_2] + param_dict["LIST_INDEX_L2_BY_NMODE2"][i] = index_l2_unique + param_dict["LIST_NMODE2_BY_INDEX_L2"][index_l2_unique].append(i) + except: + print("WARNING: the list of noise 2 L-operators contains an L-operator " + "not associated with any existing thermal environment. The noise " + "associated with this L-operator will be discarded!") + + return param_dict \ No newline at end of file diff --git a/src/mesohops/eom/eom_functions.py b/src/mesohops/eom/eom_functions.py index f53245c..9ab6e52 100644 --- a/src/mesohops/eom/eom_functions.py +++ b/src/mesohops/eom/eom_functions.py @@ -214,15 +214,16 @@ def calc_norm_corr( def calc_LT_corr( list_LT_coeff, list_L2, list_avg_L2, list_L2_sq): - """ + r""" Computes the low-temperature correction factor associated with each member of the hierarchy in the nonlinear equation of motion. The factor is given by the sum over - the low-temperature correction coefficients and associated L-operators c_n and L_n: + the low-temperature correction coefficients and associated L-operators ``c_n`` and + ``L_n``: \sum_n conj(c_n)L_n to all auxiliary wave functions and \sum_n c_n( - L_n)L_n - to the physical wave function, where c_n is the nth low-temperature correction - factor, and L_n is the nth L-operator associated with that factor. + to the physical wave function, where ``c_n`` is the nth low-temperature correction + factor, and ``L_n`` is the nth L-operator associated with that factor. Parameters ---------- @@ -257,13 +258,13 @@ def calc_LT_corr( def calc_LT_corr_to_norm_corr( list_LT_coeff, list_avg_L2, list_avg_L2_sq ): - """ + r""" Computes the low-temperature correction to the normalization factor in the normalized nonlinear equation of motion. The correction is given by the sum over - the low-temperature correction coefficients and associated L-operators c_n and L_n: + the low-temperature correction coefficients and associated L-operators ``c_n`` and ``L_n``: \sum_n Re[c_n](2^2 - ), - where c_n is the nth low-temperature correction - factor, and L_n is the nth L-operator associated with that factor. + where ``c_n`` is the nth low-temperature correction + factor, and ``L_n`` is the nth L-operator associated with that factor. Parameters ---------- @@ -285,12 +286,12 @@ def calc_LT_corr_to_norm_corr( def calc_LT_corr_linear( list_LT_coeff, list_L2_sq ): - """ + r""" Computes the low-temperature correction factor associated with each member of the hierarchy in the linear equation of motion. The factor is given by the sum over the - low-temperature correction coefficients and associated L-operators c_n and L_n: + low-temperature correction coefficients and associated L-operators ``c_n`` and ``L_n``: -\sum_n c_nL_n^2, - where c_n is the nth low-temperature correction factor, and L_n is + where ``c_n`` is the nth low-temperature correction factor, and ``L_n`` is the nth L-operator associated with that factor. NOTE: This correction should only be applied to the physical wavefunction. diff --git a/src/mesohops/eom/hops_eom.py b/src/mesohops/eom/hops_eom.py index 01ca2e3..aa0487d 100644 --- a/src/mesohops/eom/hops_eom.py +++ b/src/mesohops/eom/hops_eom.py @@ -14,8 +14,8 @@ from mesohops.util.exceptions import UnsupportedRequest __title__ = "Equations of Motion" -__author__ = "D. I. G. Bennett, B. Citty" -__version__ = "1.2" +__author__ = "D. I. G. B. Raccah, B. Citty" +__version__ = "1.6" EOM_DICT_DEFAULT = { "TIME_DEPENDENCE": False, @@ -30,11 +30,13 @@ EOM_DICT_TYPES = { "TIME_DEPENDENCE": [type(False)], - "EQUATION_OF_MOTION": [type(str())], + "EQUATION_OF_MOTION": [str], "ADAPTIVE": [type(False)], - "DELTA_A": [type(1.0)], - "DELTA_S": [type(1.0)], - "UPDATE_STEP": [type(1.0), type(False)], + "ADAPTIVE_H": [type(False)], + "ADAPTIVE_S": [type(False)], + "DELTA_A": [type(1.0), type(1)], + "DELTA_S": [type(1.0), type(1)], + "UPDATE_STEP": [type(1.0), type(False), None, type(1)], "F_DISCARD": [type(0.0)] } @@ -412,7 +414,10 @@ def dsystem_dt( Φ_view_red = Φ_view_F[list_L2_masks[rel_index][1],:] list_L2_csr_red = list_L2_csr[rel_index][list_L2_masks[rel_index][2]] - Φ_deriv_view_F[list_L2_masks[rel_index][0],:] += (z_hat1_tmp[j] - 2.0j * np.real(z_tmp2[j])) * (list_L2_csr_red @ Φ_view_red) + Φ_deriv_view_F[list_L2_masks[rel_index][0],:] += ( + (z_hat1_tmp[j] - 1.0j * z_tmp2[j]) * + (list_L2_csr_red @ Φ_view_red) + ) Φ_view_red = Φ_view_C[list_hier_mask_Zp1[rel_index][1],:] Z2_kp1_red = Z2_kp1[rel_index][list_hier_mask_Zp1[rel_index][2]] @@ -527,7 +532,7 @@ def dsystem_dt( Φ_deriv_view = np.asarray(Φ_deriv).reshape([system.size,hierarchy.size],order="F") Φ_view = np.asarray(Φ).reshape([system.size,hierarchy.size],order="F") for j in range(len(list_L2_csr)): - Φ_deriv_view += (z_hat1_tmp[j] - 2.0j * np.real(z_rnd2_tmp[j])) * ( + Φ_deriv_view += (z_hat1_tmp[j] - 1.0j * z_rnd2_tmp[j]) * ( list_L2_csr[j] @ Φ_view) # Calculates dz/dt diff --git a/src/mesohops/integrator/integrator_rk.py b/src/mesohops/integrator/integrator_rk.py index 813ec07..40f0e7a 100644 --- a/src/mesohops/integrator/integrator_rk.py +++ b/src/mesohops/integrator/integrator_rk.py @@ -26,8 +26,9 @@ def runge_kutta_step(dsystem_dt, phi, z_mem, z_rnd, z_rnd2, tau): Random numbers for the bath (at three time points) [units: cm^-1]. 5. z_rnd2 : np.array(complex) - Secondary real contribution to the noise (at three time points). - Imaginary portion discarded in dsystem_dt [units: cm^-1]. + Secondary, (typically) real contribution to the noise (at three time + points). Imaginary portion discarded by the FLAG_REAL key of the + noise object's parameter dictionary [units: cm^-1]. For primary use-case, see: "Exact open quantum system dynamics using the Hierarchy of Pure States diff --git a/src/mesohops/noise/hops_noise.py b/src/mesohops/noise/hops_noise.py index 75c3257..f29e2cd 100644 --- a/src/mesohops/noise/hops_noise.py +++ b/src/mesohops/noise/hops_noise.py @@ -10,8 +10,8 @@ from mesohops.util.physical_constants import precision # constant __title__ = "Pyhops Noise" -__author__ = "D. I. G. Bennett, J. K. Lynd" -__version__ = "1.2" +__author__ = "D. I. G. B. Raccah, B. Citty, J. K. Lynd" +__version__ = "1.6" # NOISE MODELS: # ============= @@ -26,7 +26,8 @@ "RAND_MODEL": "SUM_GAUSSIAN", # SUM_GAUSSIAN or BOX_MULLER "STORE_RAW_NOISE": False, "NOISE_WINDOW": None, - "ADAPTIVE": False + "ADAPTIVE": False, + "FLAG_REAL": False } NOISE_TYPE_DEFAULT = { @@ -38,7 +39,8 @@ "RAND_MODEL": [str], "STORE_RAW_NOISE": [type(False)], "NOISE_WINDOW": [type(None), type(1.0), type(1)], - "ADAPTIVE": [bool] + "ADAPTIVE": [bool], + "FLAG_REAL": [bool], } @@ -216,14 +218,7 @@ def _prepare_noise(self, new_lop): if self.param['STORE_RAW_NOISE']: print("Raw noise is identical to correlated noise in the ZERO noise " "model.") - z_correlated = np.zeros([n_l2, n_taus], dtype=np.complex64) - if self.param['INTERPOLATE']: - self._noise = interp1d(self.param['T_AXIS'], z_correlated, kind='cubic', - axis=1) - elif self.param['ADAPTIVE']: - new_noise = np.complex64(z_correlated) - else: - self._noise = np.complex64(z_correlated) + self._noise = 0 # FFTfilter case: elif self.param["MODEL"] == "FFT_FILTER": @@ -318,13 +313,17 @@ def _prepare_noise(self, new_lop): #Add new noise to self._noise if self.param['ADAPTIVE']: - for (i,lop) in enumerate(new_lop): - self._row += [lop]*n_taus - self._col += list(np.arange(n_taus)) - self._data += list(new_noise[i,:]) - self._noise = sp.sparse.coo_array((self._data,(self._row,self._col)), - shape=(self.param['N_L2'],len(self.param["T_AXIS"])), - dtype=np.complex64).tocsc() + # Leave the noise as a 0 integer for noise model ZERO. + if self.param['MODEL'] == 'ZERO': + pass + else: + for (i,lop) in enumerate(new_lop): + self._row += [lop]*n_taus + self._col += list(np.arange(n_taus)) + self._data += list(new_noise[i,:]) + self._noise = sp.sparse.coo_array((self._data,(self._row,self._col)), + shape=(self.param['N_L2'],len(self.param["T_AXIS"])), + dtype=np.complex64).tocsc() # Update lop_active so get_noise knows when to call prepare_noise self._lop_active = list(set(self._lop_active) | set(new_lop)) @@ -354,6 +353,12 @@ def get_noise(self, t_axis, list_lop=None): 1. Z2_noise : np.array 2D array of noise values, shape (list_lop, t_axis) sampled at the given time points. """ + if list_lop is None: + list_lop = np.arange(self.param["N_L2"]) + + if self.param["MODEL"] == "ZERO": + return np.zeros([len(list_lop), len(t_axis)], dtype=np.complex64) + if self._noise is None: if self.param["ADAPTIVE"]: self._noise = sp.sparse.coo_array( (self.param['N_L2'], @@ -363,8 +368,6 @@ def get_noise(self, t_axis, list_lop=None): self._noise = np.zeros([self.param["N_L2"], len(self.param["T_AXIS"])], dtype=np.complex64) - if list_lop is None: - list_lop = np.arange(self.param["N_L2"]) new_lop = list(set(list_lop) - set(self._lop_active)) #Prepare noise for all L-operators not already prepared. if len(new_lop) > 0: @@ -377,6 +380,8 @@ def get_noise(self, t_axis, list_lop=None): if self.param["NOISE_WINDOW"] is not None: print("Warning: noise windowing is not supported while using " "interpolated noise.") + if self.param["FLAG_REAL"]: + return np.real(self._noise(t_axis)) return self._noise(t_axis) else: @@ -440,7 +445,9 @@ def get_noise(self, t_axis, list_lop=None): "Off axis t-samples when INTERPOLATE = False", "NoiseModel.get_noise()", ) - + if self.param["FLAG_REAL"]: + return np.real(self._noise_to_array(self.Z2_windowed, it_list, + list_lop)) return self._noise_to_array(self.Z2_windowed, it_list, list_lop) def _prepare_rand(self,new_lop=None): diff --git a/src/mesohops/noise/prepare_functions.py b/src/mesohops/noise/prepare_functions.py index f7a6a8e..729704d 100644 --- a/src/mesohops/noise/prepare_functions.py +++ b/src/mesohops/noise/prepare_functions.py @@ -1,12 +1,13 @@ +import warnings from collections import Counter from mesohops.noise.hops_noise import HopsNoise from mesohops.util.helper_functions import array_to_tuple __title__ = "preparation functions" -__author__ = "D. I. G. Bennett, J. K. Lynd" -__version__ = "1.2" +__author__ = "D. I. G. B. Raccah, J. K. Lynd" +__version__ = "1.6" -def prepare_hops_noise(noise_param, system_param, flag=1): +def prepare_hops_noise(noise_param, system_param, noise_type=1): """ Returns the proper noise class given the user inputs. @@ -18,8 +19,11 @@ def prepare_hops_noise(noise_param, system_param, flag=1): 2. system_param : dict Dictionary of system parameters. - 3. flag : int - 1-NOISE1 parameters, 2-NOISE2 parameters. + 3. noise_type : int + Defines whether the noise is noise 1 (the noise associated with + the bath correlation function modes of the hierarchy) or noise 2 + (an explicit time-dependence in the system Hamiltonian), which is + typically purely real (options: 1 or 2). Returns ------- @@ -27,7 +31,8 @@ def prepare_hops_noise(noise_param, system_param, flag=1): """ # DETERMINE CORRELATION PARAMETERS - if flag == 1: + # Skips checking for all the noise 2 params if noise 2 is defaulted to a zero model. + if noise_type == 1 or noise_param["MODEL"] == "ZERO": # Creates list of unique l2 tuples in order they appear in "L_NOISE1" l2_as_tuples = [array_to_tuple(L2) for L2 in system_param["L_NOISE1"]] list_unique_l2_as_tuples = list(Counter(l2_as_tuples)) @@ -38,7 +43,17 @@ def prepare_hops_noise(noise_param, system_param, flag=1): "NMODE_BY_LIND": system_param["LIST_NMODE1_BY_INDEX_L2"], "CORR_PARAM": system_param["PARAM_NOISE1"], } - elif flag == 2: + + if not "FLAG_REAL" in noise_param.keys(): + noise_param["FLAG_REAL"] = False + if noise_param["FLAG_REAL"] == True: + warnings.warn("Noise 1 should never be flagged real. For a purely " + "real noise, set PARAM_NOISE1 to ensure that noise 1 " + "goes to zero and use noise 2 for the real noise " + "instead.") + noise_param["FLAG_REAL"] = False + + elif noise_type == 2: # Creates list of unique l2 tuples in order they appear in "L_NOISE2" l2_as_tuples = [array_to_tuple(L2) for L2 in system_param["L_NOISE2"]] list_unique_l2_as_tuples = list(Counter(l2_as_tuples)) @@ -49,8 +64,20 @@ def prepare_hops_noise(noise_param, system_param, flag=1): "NMODE_BY_LIND": system_param["LIST_NMODE2_BY_INDEX_L2"], "CORR_PARAM": system_param["PARAM_NOISE2"], } + # Noise 2 defaults to purely-real noise, in accordance with the most common + # use cases: time-dependence of the Hermitian system Hamiltonian and the + # noise 2 laid out in "Exact open quantum system dynamics using the Hierarchy + # of Pure States (HOPS)." Richard Hartmann and Walter T. Strunz J. Chem. + # Theory Comput. 13, p. 5834-5845 (2017). + if "FLAG_REAL" not in noise_param.keys(): + noise_corr["FLAG_REAL"] = True + warnings.warn("Noise 2 FLAG_REAL not specified: setting to True. If noise 2 " + "should not be purely real, specify FLAG_REAL as False in " + "the noise 2 parameters during initialization of the HOPS " + "trajectory.") + else: - Exception("Unknown flag value in prepare_hops_noise") + Exception("Unknown noise_type value in prepare_hops_noise") # Instantiate a HopsNoise subclass # -------------------------------- diff --git a/src/mesohops/storage/hops_storage.py b/src/mesohops/storage/hops_storage.py index c65e9d1..2d5c1e4 100644 --- a/src/mesohops/storage/hops_storage.py +++ b/src/mesohops/storage/hops_storage.py @@ -3,10 +3,11 @@ from mesohops.storage.storage_functions import storage_default_func as default_func from mesohops.util.exceptions import UnsupportedRequest +from mesohops.util.git_utils import get_git_commit_hash __title__ = "Storage Class" -__author__ = "D. I. G. Bennett, L. Varvelo" -__version__ = "1.2" +__author__ = "D. I. G. Bennett, L. Varvelo, J. K. Lynd" +__version__ = "1.6" class HopsStorage: @@ -16,14 +17,17 @@ class HopsStorage: __slots__ = ( # --- Core basis components --- - '_adaptive', # Adaptive calculation flag - '_n_dim', # System dimension + '_adaptive', # Adaptive calculation flag + '_n_dim', # System dimension # --- Storage containers --- - 'storage_dic', # Main storage dictionary (current data) - 'dic_save', # Save function dictionary + 'storage_dic', # Main storage dictionary (current data) + 'dic_save', # Save function dictionary 'data', # Data storage 'metadata', # Metadata dictionary + + # --- Storage managers --- + 'storage_time', # Controls which time points are saved ) def __init__(self, adaptive, storage_dic): @@ -40,12 +44,22 @@ def __init__(self, adaptive, storage_dic): self._adaptive = False self._n_dim = 0 self.storage_dic = storage_dic + if "STORAGE_TIME" in storage_dic: + self.storage_time = storage_dic["STORAGE_TIME"] + del storage_dic["STORAGE_TIME"] + else: + self.storage_time = True self.dic_save = {} self.data = {} self.adaptive = adaptive - self.metadata = {"INITIALIZATION_TIME": 0, - "LIST_PROPAGATION_TIME": []} - + + # Initialize metadata with git commit hash + git_commit = get_git_commit_hash() + self.metadata = { + "INITIALIZATION_TIME": 0, + "LIST_PROPAGATION_TIME": [], + "GIT_COMMIT_HASH": git_commit + } def __repr__(self): key_dict = [] @@ -147,6 +161,45 @@ def store_step(self, **kwargs): for (key,value) in self.dic_save.items(): self.data[key].append(self.dic_save[key](**kwargs)) + def check_storage_time(self, t): + """ + Checks whether the current time is a time point that the user indicated + should be saved. If the storage time parameter is a list or array, returns + true if the current time is an element of the iterable. If the storage time is + a number, returns true if the current time is an integer multiple of that + number. Finally, a boolean storage time indicates that either all (True) or no + (False) time points should be saved. + + Parameters + ---------- + 1. t : float + The current time [units: fs]. + + Returns + ------- + 1. save_flag : bool + Whether the current time point's data should be saved. + + """ + if isinstance(self.storage_time, (list, np.ndarray)): + return t in self.storage_time + + elif (isinstance(self.storage_time, (int, float)) and not + isinstance(self.storage_time, bool)): + return (t%self.storage_time == 0) + + elif isinstance(self.storage_time, bool): + return self.storage_time + + else: + raise UnsupportedRequest("Type {} is not allowed in the " + "STORAGE_TIME key of storage parameters " + "for checking whether the current time " + "should be saved". + format(type(self.storage_time)), + 'check_storage_time', + override=True) + @property def n_dim(self): return self._n_dim @@ -154,4 +207,3 @@ def n_dim(self): @n_dim.setter def n_dim(self, N_dim): self._n_dim = N_dim - diff --git a/src/mesohops/timing/timing_models/absorption.py b/src/mesohops/timing/timing_models/absorption.py index c119072..370d3c0 100644 --- a/src/mesohops/timing/timing_models/absorption.py +++ b/src/mesohops/timing/timing_models/absorption.py @@ -6,7 +6,7 @@ from mesohops.trajectory.dyadic_spectra import (prepare_spectroscopy_input_dict, prepare_chromophore_input_dict, prepare_convergence_parameter_dict) -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_spectroscopy_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -41,7 +41,7 @@ M2_mu_ge = np.tile(np.array([0, 0, 1]), (nsite, 1)) H2_sys_hamiltonian = generate_spectroscopy_hamiltonian(nsite, V) list_loperators = generate_spectroscopy_loperators(nsite) -list_modes = bcf_convert_sdl_to_exp(e_lambda, gamma, 0, temp) +list_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp) chromophore_input = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict={"list_lop": list_loperators, diff --git a/src/mesohops/timing/timing_models/fluorescence.py b/src/mesohops/timing/timing_models/fluorescence.py index 8a39239..7fa98ab 100644 --- a/src/mesohops/timing/timing_models/fluorescence.py +++ b/src/mesohops/timing/timing_models/fluorescence.py @@ -6,7 +6,7 @@ from mesohops.trajectory.dyadic_spectra import (prepare_spectroscopy_input_dict, prepare_chromophore_input_dict, prepare_convergence_parameter_dict) -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_spectroscopy_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -42,7 +42,7 @@ M2_mu_ge = np.tile(np.array([0, 0, 1]), (nsite, 1)) H2_sys_hamiltonian = generate_spectroscopy_hamiltonian(nsite, V) list_loperators = generate_spectroscopy_loperators(nsite) -list_modes = bcf_convert_sdl_to_exp(50, 50, 0, 300) +list_modes = bcf_convert_dl_to_exp(50, 50, 300) chromophore_input = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict={"list_lop": list_loperators, diff --git a/src/mesohops/timing/timing_models/holstein_1_particle.py b/src/mesohops/timing/timing_models/holstein_1_particle.py index 4c611eb..f19ea4f 100644 --- a/src/mesohops/timing/timing_models/holstein_1_particle.py +++ b/src/mesohops/timing/timing_models/holstein_1_particle.py @@ -4,7 +4,7 @@ import numpy as np from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_1_particle_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -36,7 +36,7 @@ list_loperators = generate_holstein_1_particle_loperators(nsite) # Gets the list of bath correlation function modes for each independent environment -list_dl_modes = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) +list_dl_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp) def prepare_gw_sysbath(list_lop, list_modes): """ diff --git a/src/mesohops/timing/timing_models/holstein_2_particle.py b/src/mesohops/timing/timing_models/holstein_2_particle.py index eb7f4a9..844065b 100644 --- a/src/mesohops/timing/timing_models/holstein_2_particle.py +++ b/src/mesohops/timing/timing_models/holstein_2_particle.py @@ -4,7 +4,7 @@ import numpy as np from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_2_particle_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -39,7 +39,7 @@ list_loperators = generate_holstein_2_particle_loperators(nsite) # Gets the list of bath correlation function modes for each independent environment -list_dl_modes = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) +list_dl_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp) def prepare_gw_sysbath(list_lop, list_modes): """ diff --git a/src/mesohops/timing/timing_models/longedge_filter.py b/src/mesohops/timing/timing_models/longedge_filter.py index a30d0a1..1f008e7 100644 --- a/src/mesohops/timing/timing_models/longedge_filter.py +++ b/src/mesohops/timing/timing_models/longedge_filter.py @@ -4,7 +4,7 @@ import numpy as np from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp_with_Matsubara +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_1_particle_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -38,7 +38,7 @@ list_loperators = generate_holstein_1_particle_loperators(nsite) # Gets the list of bath correlation function modes for each independent environment -list_dl_modes = bcf_convert_dl_to_exp_with_Matsubara(e_lambda, gamma, temp, k_mats) +list_dl_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp, k_mats) def prepare_gw_sysbath(list_lop, list_modes): diff --git a/src/mesohops/timing/timing_models/markovian_filter.py b/src/mesohops/timing/timing_models/markovian_filter.py index dcc8a30..33c9379 100644 --- a/src/mesohops/timing/timing_models/markovian_filter.py +++ b/src/mesohops/timing/timing_models/markovian_filter.py @@ -4,7 +4,7 @@ import numpy as np from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp_with_Matsubara +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_1_particle_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -37,7 +37,7 @@ list_loperators = generate_holstein_1_particle_loperators(nsite) # Gets the list of bath correlation function modes for each independent environment -list_dl_modes = bcf_convert_dl_to_exp_with_Matsubara(e_lambda, gamma, temp, k_mats) +list_dl_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp, k_mats) def prepare_gw_sysbath(list_lop, list_modes): diff --git a/src/mesohops/timing/timing_models/peierls.py b/src/mesohops/timing/timing_models/peierls.py index da0d558..2174b38 100644 --- a/src/mesohops/timing/timing_models/peierls.py +++ b/src/mesohops/timing/timing_models/peierls.py @@ -4,7 +4,7 @@ import numpy as np from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_1_particle_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -36,7 +36,7 @@ list_loperators = generate_peierls_1_particle_loperators(nsite) # Gets the list of bath correlation function modes for each independent environment -list_dl_modes = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) +list_dl_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp) def prepare_gw_sysbath(list_lop, list_modes): """ diff --git a/src/mesohops/timing/timing_models/triangular_filter.py b/src/mesohops/timing/timing_models/triangular_filter.py index eeb2b79..9191f12 100644 --- a/src/mesohops/timing/timing_models/triangular_filter.py +++ b/src/mesohops/timing/timing_models/triangular_filter.py @@ -4,7 +4,7 @@ import numpy as np from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp_with_Matsubara +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.timing.helper_functions.hamiltonian_generation import ( generate_1_particle_hamiltonian) from mesohops.timing.helper_functions.loperator_generation import ( @@ -37,7 +37,7 @@ list_loperators = generate_holstein_1_particle_loperators(nsite) # Gets the list of bath correlation function modes for each independent environment -list_dl_modes = bcf_convert_dl_to_exp_with_Matsubara(e_lambda, gamma, temp, k_mats) +list_dl_modes = bcf_convert_dl_to_exp(e_lambda, gamma, temp, k_mats) def prepare_gw_sysbath(list_lop, list_modes): diff --git a/src/mesohops/trajectory/dyadic_spectra.py b/src/mesohops/trajectory/dyadic_spectra.py index 664e63e..a516955 100644 --- a/src/mesohops/trajectory/dyadic_spectra.py +++ b/src/mesohops/trajectory/dyadic_spectra.py @@ -8,6 +8,7 @@ __author__ = "D. I. G. B. Raccah, A. Hartzell, T. Gera, J. K. Lynd" __version__ = "1.5" + class DyadicSpectra(DyadicTrajectory): """ Acts as an interface to calculate spectra using the Dyadic HOPS method. @@ -15,47 +16,47 @@ class DyadicSpectra(DyadicTrajectory): __slots__ = ( # --- Initialization and tracking --- - '__initialized', # Initialization status flag + '__initialized', # Initialization status flag # --- Spectroscopy parameters --- - 'spectrum_type', # Type of spectrum to calculate - 't_1', # First propagation time - 't_2', # Second propagation time - 't_3', # Third propagation time - 'list_t', # List of propagation times - 'E_1', # First field definition - 'E_2', # Second field definition - 'E_3', # Third field definition - 'E_sig', # Signal field definition - 'list_ket_sites', # Ket sites excited by field - 'list_bra_sites', # Bra sites excited by field + 'spectrum_type', # Type of spectrum to calculate + 't_1', # First propagation time + 't_2', # Second propagation time + 't_3', # Third propagation time + 'list_t', # List of propagation times + 'E_1', # First field definition + 'E_2', # Second field definition + 'E_3', # Third field definition + 'E_sig', # Signal field definition + 'list_ket_sites', # Ket sites excited by field + 'list_bra_sites', # Bra sites excited by field # --- Chromophore parameters --- - 'M2_mu_ge', # Transition dipole matrix - 'n_chromophore', # Number of chromophores + 'M2_mu_ge', # Transition dipole matrix + 'n_chromophore', # Number of chromophores 'H2_sys_hamiltonian', # System Hamiltonian - 'lop_list_hier', # L-operators associated with hierarchy modes - 'gw_sysbath_hier', # Hierarchy mode parameters - 'lop_list_noise', # L-operators associated with noise - 'gw_sysbath_noise', # Noise mode parameters - 'lop_list_ltc', # L-operators associated with LTC - 'ltc_param', # Low-temperature correction parameters + 'lop_list_hier', # L-operators associated with hierarchy modes + 'gw_sysbath_hier', # Hierarchy mode parameters + 'lop_list_noise', # L-operators associated with noise + 'gw_sysbath_noise', # Noise mode parameters + 'lop_list_ltc', # L-operators associated with LTC + 'ltc_param', # Low-temperature correction parameters # --- Convergence parameters --- - 't_step', # Time step - 'max_hier', # Maximum hierarchy depth - 'delta_a', # Auxiliary derivative error bound - 'delta_s', # State derivative error bound - 'set_update_step', # Update step - 'set_f_discard', # Discard fraction + 't_step', # Time step + 'max_hier', # Maximum hierarchy depth + 'delta_a', # Auxiliary derivative error bound + 'delta_s', # State derivative error bound + 'set_update_step', # Update step + 'set_f_discard', # Discard fraction 'static_filter_list', # Static hierarchy filters # --- State dimensions --- - 'n_state_hilb', # Hilbert space dimension - 'n_state_dyad', # Dyadic space dimension + 'n_state_hilb', # Hilbert space dimension + 'n_state_dyad', # Dyadic space dimension # --- Noise configuration --- - 'noise_param' # Noise parameters + 'noise_param' # Noise parameters ) def __init__(self, spectroscopy_dict, chromophore_dict, convergence_dict, seed): @@ -142,6 +143,9 @@ def __init__(self, spectroscopy_dict, chromophore_dict, convergence_dict, seed): "TLEN": float(1000 + np.sum(self.list_t)), "TAU": 0.5 if self.t_step % 1 == 0 else self.t_step / 2} + # Noise 2 not currently supported. + self.noise2_param = None + # Preparing hierarchy parameter dictionary hierarchy_param = {"MAXHIER": self.max_hier} if self.static_filter_list: @@ -151,8 +155,8 @@ def __init__(self, spectroscopy_dict, chromophore_dict, convergence_dict, seed): storage_param = {} # Initializing DyadicTrajectory class - super().__init__(system_param, eom_param, self.noise_param, hierarchy_param, - storage_param) + super().__init__(system_param, eom_param, self.noise_param, + self.noise2_param, hierarchy_param, storage_param) def initialize(self): """ @@ -560,20 +564,32 @@ def prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict): Brian Citty, Jacob K. Lynd, et al. J. Chem. Phys. 160, 144118 (2024) - e. static_filter_list: list, optional - List of static filters applied to the - hierarchy. Each filter is defined by a list - of the form [filter_name, filter_params]. - OPTIONS: - -------- - 1. "Markovian": auxiliary wave functions + e. static_filter_list: list(list), optional + List of static filters applied to the + hierarchy. Each filter is defined by a list + of the form [filter_name, filter_params]. + This means the full structure of + static_filter_list is: + [filter1, filter2, ...] where each + filter is a list of the form + [filter_name, filter_params] + where filter_name is a string and + filter_params is a list of parameters + defining the filter. The length of the + boolean list in filter_params should match + the number of modes in list_modes or + list_modes_by_bath, depending on which + is defined in bath_dict. + OPTIONS: + -------- + 1. "Markovian": auxiliary wave functions associated with filtered modes are only included in the hierarchy if they are depth 1. filter_params should be a list of booleans: True for filtered modes, False for unfiltered modes. - 2. "Triangular": auxiliary wave functions + 2. "Triangular": auxiliary wave functions associated with filtered modes are only included in the hierarchy if they are at or below depth kmax2. filter_params = @@ -583,7 +599,7 @@ def prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict): False for unfiltered modes. Note kmax2 is an integer. - 3. "LongEdge": auxiliary wave functions + 3. "LongEdge": auxiliary wave functions associated with filtered modes are only included in the hierarchy if they are at or below depth kmax2 OR only have depth @@ -594,23 +610,18 @@ def prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict): False for unfiltered modes. Note kmax2 is an integer. - If the list of booleans determining which - modes are filtered does not match the full - set of modes in the "GW_SYSBATH" key of the - system parameter dictionary, the trajectory - will raise an error and the simulation will - be terminated. The use of any static - hierarchy filter may reduce the accuracy of a - given calculation; that is, static hierarchy - filters must be tested like any other - convergence parameters. - - For more details on static filters, see: - "MesoHOPS: Size-invariant scaling - calculations of multi-excitation open quantum - systems." - Brian Citty, Jacob K. Lynd, et al. - J. Chem. Phys. 160, 144118 (2024) + The use of any static + hierarchy filter may reduce the accuracy of a + given calculation; that is, static hierarchy + filters must be tested like any other + convergence parameters. + + For more details on static filters, see: + "MesoHOPS: Size-invariant scaling + calculations of multi-excitation open quantum + systems." + Brian Citty, Jacob K. Lynd, et al. + J. Chem. Phys. 160, 144118 (2024) Returns ------- @@ -618,177 +629,223 @@ def prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict): Dictionary of chromophore parameters needed for DyadicSpectra class. """ - # Defining number of chromophores + # Define number of chromophores and validate M2_mu_ge structure n_chromophore = len(M2_mu_ge) - - # Checking M2_mu_ge input structure M2_mu_ge = np.array(M2_mu_ge) if M2_mu_ge.shape[1] != 3: raise ValueError( "M2_mu_ge must be a numpy array with shape (n_chromophore, 3).") - # Setting default nmodes_LTC to 0 and checking nmodes_LTC input structure + # Clean up bath_dict: convert arrays to lists and remove None/0 values + for key, value in list(bath_dict.items()): + # Convert numpy arrays to lists for consistency + if type(value) == np.ndarray: + bath_dict[key] = list(value) + # Remove keys with None/0 values (except nmodes_LTC) + elif (key != "nmodes_LTC") and (value is None or value == 0): + del bath_dict[key] + + # Set default value for nmodes_LTC if not provided or None if "nmodes_LTC" not in bath_dict.keys() or bath_dict["nmodes_LTC"] is None: bath_dict["nmodes_LTC"] = 0 - + # Validate nmodes_LTC if provided elif not isinstance(bath_dict["nmodes_LTC"], int): raise ValueError("nmodes_LTC must be an integer or None.") - elif bath_dict["nmodes_LTC"] < 0: raise ValueError("nmodes_LTC must be >= 0 or set to None.") - # Checking static_filter_list input structure - if "static_filter_list" in bath_dict.keys(): - if not isinstance(bath_dict["static_filter_list"], list): - raise ValueError("static_filter_list must be a list.") - - elif len(bath_dict["static_filter_list"]) != 2: - raise ValueError("static_filter_list must be a 2-element list of the form: " - "[filter_name, filter_params].") - - elif bath_dict["static_filter_list"][0] not in ["Markovian", "Triangular", - "LongEdge"]: - raise ValueError("Filter name must be 'Markovian', 'Triangular', or " - "'LongEdge'.") - - elif (bath_dict["static_filter_list"][0] == "Markovian" and - not all(isinstance(boolean, bool) for boolean in - bath_dict["static_filter_list"][1])): - raise ValueError( - "filter_params for Markovian filter must be a list of booleans.") - - elif (bath_dict["static_filter_list"][0] in ["Triangular", "LongEdge"] and - len(bath_dict["static_filter_list"][1]) != 2): - raise ValueError("The Triangular and LongEdge filter_params must be a " - "list of booleans and an integer.") - - elif (bath_dict["static_filter_list"][0] in ["Triangular", "LongEdge"] and - not isinstance(bath_dict["static_filter_list"][1][1], int)): - raise ValueError("The second entry in filter_params for the Triangular " - "and LongEdge filters must be an integer.") - - elif bath_dict["static_filter_list"][0] in ["Triangular", "LongEdge"]: - if not all(isinstance(boolean, bool) for boolean in - bath_dict["static_filter_list"][1][0]): - raise ValueError("The first entry in filter_params for the " - "Triangular and LongEdge filters must be a " - "list of booleans.") - - elif (bath_dict["static_filter_list"][0] in ["Triangular", "LongEdge"] and - bath_dict["static_filter_list"][1][1] < 0): - raise ValueError("Triangular and LongEdge filter_params must have a " - "positive integer as the second element.") - - # Checking bath_dict input structure - for key, value in list(bath_dict.items()): - # Setting bath_dict arrays to lists - if type(value) == np.ndarray: - bath_dict[key] = list(value) - - # Removing keys with None/0 values from bath_dict (except nmodes_LTC) - elif (key != "nmodes_LTC") and (value is None or value == 0): - del bath_dict[key] + # Check that users don't try to use LTC with static filters + if bath_dict["nmodes_LTC"] > 0 and "static_filter_list" in bath_dict: + raise ValueError("The use of static hierarchy filters with low-temperature " + "correction is not currently supported.") - # Checking list_modes/list_modes_by_bath over-definition + # Check that list_modes and list_modes_by_bath are not both defined if "list_modes_by_bath" in bath_dict.keys() and "list_modes" in bath_dict.keys(): raise ValueError( "list_modes_by_bath and list_modes should not both be defined.") - # Setting default list_lop if not defined + # Set default list_lop if not defined (site-projection operators) if "list_lop" not in bath_dict.keys(): - # Site-projection L-operators bath_dict["list_lop"] = [sparse.coo_matrix(([1], ([chrom + 1], [chrom + 1])), shape=(n_chromophore + 1, n_chromophore + 1)) for chrom in range(n_chromophore)] - # Checking list_modes_by_bath/list_modes structure - if "list_modes_by_bath" in bath_dict: - # Checking list_modes_by_bath/list_lop compatibility + # Process list_modes if provided + if "list_modes" in bath_dict: + # Validate list_modes structure (must have paired Gs and Ws) + if len(bath_dict["list_modes"]) % 2 != 0: + raise ValueError("list_modes should contain paired Gs and Ws, which " + "guarantees an even number of elements.") + + # Create list_modes_by_bath by repeating list_modes for each bath + bath_dict["list_modes_by_bath"] = [bath_dict["list_modes"] for _ in + range(len(bath_dict["list_lop"]))] + + # Process list_modes_by_bath if provided + elif "list_modes_by_bath" in bath_dict: + # Validate compatibility with list_lop if len(bath_dict["list_modes_by_bath"]) != len(bath_dict["list_lop"]): raise ValueError( "list_modes_by_bath and list_lop must have the same length.") - # Checking that list_modes_by_bath is a list of lists of paired Gs and Ws + # Validate structure of each sublist (must have paired Gs and Ws) for sublist in bath_dict["list_modes_by_bath"]: if not isinstance(sublist, list): raise ValueError("list_modes_by_bath must be a list of lists.") elif len(sublist) % 2 != 0: - raise ValueError("list_modes_by_bath should contain paired Gs and Ws, " - "which guarantees an even number of elements in each " - "sublist.") + raise ValueError("sublists within list_modes_by_bath should contain " + "paired Gs and Ws, which guarantees an even number of " + "elements in each sublist.") - elif "list_modes" in bath_dict: - # Checking that list_modes is a list of paired Gs and Ws - if len(bath_dict["list_modes"]) % 2 != 0: - raise ValueError("list_modes should contain paired Gs and Ws, which " - "guarantees an even number of elements.") + # Require either list_modes or list_modes_by_bath is defined else: raise ValueError("Either list_modes_by_bath or list_modes must be defined.") - # Preparing empty chromophore dictionary - gw_sysbath = [] - list_lop_sysbath_by_mode = [] - gw_noise = [] - list_lop_noise_by_mode = [] - list_lop_ltc = [] - list_ltc_param = [] + # Process static_filter_list if provided + if "static_filter_list" in bath_dict.keys(): + # Validate that static_filter_list is a list + if not isinstance(bath_dict["static_filter_list"], list): + raise ValueError("static_filter_list must be a list.") - # Populating chromophore dictionary + # Validate each filter in the list + for filter_idx, filter_item in enumerate(bath_dict["static_filter_list"]): + # Check each filter is a list of the form [filter_name, filter_params] + if not isinstance(filter_item, list): + raise ValueError( + f"Error in filter {filter_idx}: static_filter_list must be a list.") + elif len(filter_item) != 2: + raise ValueError( + f"Error in filter {filter_idx}: each filter in static_filter_list " + f"must be a 2-element list of the form: " + f"[filter_name, filter_params].") + # Check filter_name is a string and is one of the allowed types + elif filter_item[0] not in ["Markovian", "Triangular", "LongEdge"]: + raise ValueError( + f"Error in filter {filter_idx}: Filter names must be 'Markovian', " + f"'Triangular', or 'LongEdge'.") + + # Validate filter parameters for Markovian filters + if filter_item[0] == "Markovian": + # Markovian filter_params should be a list of booleans + if not all(isinstance(boolean, bool) for boolean in filter_item[1]): + raise ValueError( + f"Error in filter {filter_idx}: filter_params for Markovian " + f"filters must be a list of booleans.") + + # Check that the length of filter_item[1] matches the number of modes + # in each bath if given list_modes input + if "list_modes" in bath_dict.keys(): + if len(filter_item[1]) != len(bath_dict["list_modes"]) // 2: + raise ValueError(f"Error in filter {filter_idx}: The list of " + f"booleans in filter_params must have the " + f"same length as the number of modes in " + f"each bath.") + else: + # Expand filter to cover all baths + filter_item[1] = filter_item[1] * len(bath_dict["list_lop"]) + + # If list_modes_by_bath is given, check against total number of modes + elif len(filter_item[1]) != np.sum([len(mode_list) for mode_list in + bath_dict["list_modes_by_bath"]])//2: + raise ValueError(f"Error in filter {filter_idx}: The list of " + f"booleans in filter_params must have " + f"the same length as the number of modes in all " + f"baths combined.") + + # Validate filter parameters for Triangular and LongEdge filters + elif filter_item[0] in ["Triangular", "LongEdge"]: + # Triangular and LongEdge filter_params should be a list containing a + # list of booleans and an integer + if len(filter_item[1]) != 2: + raise ValueError( + f"Error in filter {filter_idx}: Triangular/LongEdge " + f"filter_params must be a list containing a list of booleans " + f"and an integer.") + elif not isinstance(filter_item[1][1], int): + raise ValueError( + f"Error in filter {filter_idx}: The second entry in " + f"filter_params for Triangular/LongEdge filters must be an " + f"integer.") + elif not all( + isinstance(boolean, bool) for boolean in filter_item[1][0]): + raise ValueError(f"Error in filter {filter_idx}: The first entry in" + f" filter_params for Triangular/LongEdge filters " + f"must be a list of booleans.") + elif filter_item[1][1] < 0: + raise ValueError(f"Error in filter {filter_idx}: " + f"Triangular/LongEdge filter_params must have a " + f"positive integer as the second element.") + + # Check that the length of filter_item[1][0] matches the number of modes + # in each bath if given list_modes input + if "list_modes" in bath_dict.keys(): + if len(filter_item[1][0]) != len(bath_dict["list_modes"]) // 2: + raise ValueError(f"Error in filter {filter_idx}: The list of " + f"booleans in filter_params must have the same" + f" length as the number of modes in each bath.") + else: + # Expand filter to cover all baths + filter_item[1][0] = (filter_item[1][0] * + len(bath_dict["list_lop"])) + + # If list_modes_by_bath is given, check against total number of modes + elif len(filter_item[1][0]) != np.sum([len(mode_list) for mode_list in + bath_dict["list_modes_by_bath"]])//2: + raise ValueError(f"Error in filter {filter_idx}: The list of " + f"booleans in filter_params must have the same " + f"length as the number of modes in all baths " + f"combined.") + + # Initialize empty lists for chromophore dictionary components + gw_sysbath = [] # G-W tuples for hierarchy + list_lop_sysbath_by_mode = [] # L-operators for hierarchy + gw_noise = [] # G-W tuples for noise + list_lop_noise_by_mode = [] # L-operators for noise + list_lop_ltc = [] # L-operators for low-temperature correction + list_ltc_param = [] # Low-temperature correction parameters + + # Process each bath for bath in range(len(bath_dict["list_lop"])): - # Unpacking list of Gs and Ws for each bath in the list_modes_by_bath case - if "list_modes_by_bath" in bath_dict.keys(): - bath_dict["list_modes"] = bath_dict["list_modes_by_bath"][bath] - - # Converting list of Gs and Ws to list of coupled G-W tuples - list_modes_as_tuples = [(bath_dict["list_modes"][i], - bath_dict["list_modes"][i + 1]) for i in - range(0, len(bath_dict["list_modes"]), 2)] + # Convert list of Gs and Ws to list of coupled G-W tuples + list_modes_as_tuples = [(bath_dict["list_modes_by_bath"][bath][i], + bath_dict["list_modes_by_bath"][bath][i + 1]) for i in + range(0, len(bath_dict["list_modes_by_bath"][bath]), 2)] - # Checking nmodes_LTC and static_filter_list compatibility with number of modes + # Validate nmodes_LTC compatibility with number of modes if bath_dict["nmodes_LTC"] >= len(list_modes_as_tuples): - raise ValueError("nmodes_LTC must be less than the number of " - "modes in each bath.") - - elif "static_filter_list" in bath_dict.keys(): - if bath_dict["static_filter_list"][0] == 'Markovian': - if len(bath_dict["static_filter_list"][1]) != len(list_modes_as_tuples): - raise ValueError("The list of booleans in static_filter_list must " - "have the same length as the number of modes.") - - elif bath_dict["static_filter_list"][0] in ['Triangular', 'LongEdge']: - if len(bath_dict["static_filter_list"][1][0]) != len( - list_modes_as_tuples): - raise ValueError("The list of booleans in static_filter_list must " - "have the same length as the number of modes.") + raise ValueError("nmodes_LTC must be less than the number of modes in each " + "bath.") + # Initialize LTC parameter for this bath ltc_param = 0 list_lop_ltc.append(bath_dict["list_lop"][bath]) - # Appending gw_sysbath and lop_list for each G-W tuple + + # Process each mode in the bath for nmode, mode in enumerate(list_modes_as_tuples): + # Append G-W tuples and L-operators to the appropriate lists gw_noise.append(mode) list_lop_noise_by_mode.append(bath_dict["list_lop"][bath]) - # Checking if mode belongs to the low-temperature corrected modes + # Determine if mode is treated in hierarchy or with LTC if len(list_modes_as_tuples) - nmode > bath_dict["nmodes_LTC"]: - # Appending each G-W tuple treated in the hierarchy + # Mode is treated in hierarchy gw_sysbath.append(mode) list_lop_sysbath_by_mode.append(bath_dict["list_lop"][bath]) - - # Updating ltc parameter for each low-temperature corrected mode else: + # Mode is treated with low-temperature correction ltc_param += mode[0] / mode[1] - # Appending ltc parameter to the bath + # Add LTC parameter for this bath list_ltc_param.append(ltc_param) - # Returning chromophore dictionary return {"M2_mu_ge": M2_mu_ge, "n_chromophore": n_chromophore, "H2_sys_hamiltonian": H2_sys_hamiltonian, "lop_list_hier": list_lop_sysbath_by_mode, "gw_sysbath_hier": gw_sysbath, "lop_list_noise": list_lop_noise_by_mode, "gw_sysbath_noise": gw_noise, "lop_list_ltc": list_lop_ltc, "ltc_param": list_ltc_param, - "static_filter_list": bath_dict.get("static_filter_list", None)} + "static_filter_list": bath_dict.get("static_filter_list", None), + } def prepare_convergence_parameter_dict(t_step, max_hier, delta_a=0, delta_s=0, @@ -826,4 +883,3 @@ def prepare_convergence_parameter_dict(t_step, max_hier, delta_a=0, delta_s=0, return {"t_step": t_step, "max_hier": max_hier, "delta_a": delta_a, "delta_s": delta_s, "set_update_step": set_update_step, "set_f_discard": set_f_discard} - diff --git a/src/mesohops/trajectory/exp_noise.py b/src/mesohops/trajectory/exp_noise.py index c21c587..8723de0 100644 --- a/src/mesohops/trajectory/exp_noise.py +++ b/src/mesohops/trajectory/exp_noise.py @@ -3,7 +3,7 @@ def bcf_exp(t_axis, g, w): - """ + r""" This is the form of the correlation function: alpha(t) = \displastyle\ g exp(-w*t) diff --git a/src/mesohops/trajectory/hops_dyadic.py b/src/mesohops/trajectory/hops_dyadic.py index 06f9c0b..9feeea8 100644 --- a/src/mesohops/trajectory/hops_dyadic.py +++ b/src/mesohops/trajectory/hops_dyadic.py @@ -21,7 +21,8 @@ class DyadicTrajectory(HopsTrajectory): '__list_response_norm_sq', # List of normalization factors ) - def __init__(self, system_param, eom_param=None, noise_param=None, hierarchy_param=None, + def __init__(self, system_param, eom_param=None, noise_param=None, + noise2_param = None, hierarchy_param=None, storage_param=None, integration_param=None): """ Inputs @@ -76,7 +77,7 @@ def __init__(self, system_param, eom_param=None, noise_param=None, hierarchy_par self._M2_dyad_conversion(system_param[param][i])) system_param.update({param: dyad_lop_list}) - super().__init__(system_param, eom_param, noise_param, + super().__init__(system_param, eom_param, noise_param, noise2_param, hierarchy_param, storage_param, integration_param) def initialize(self, psi_ket, psi_bra, timer_checkpoint=None): @@ -135,7 +136,7 @@ def _M2_dyad_conversion(self, M2_hilbert): shape=(2*M2_dim, 2*M2_dim)) else: - M2_dyad = np.zeros((2 * M2_dim, 2 * M2_dim)) + M2_dyad = np.zeros((2 * M2_dim, 2 * M2_dim),dtype=np.complex128) M2_dyad[M2_dim:, M2_dim:] = M2_hilbert M2_dyad[:M2_dim, :M2_dim] = M2_hilbert diff --git a/src/mesohops/trajectory/hops_trajectory.py b/src/mesohops/trajectory/hops_trajectory.py index 96bffd9..2b90e0c 100644 --- a/src/mesohops/trajectory/hops_trajectory.py +++ b/src/mesohops/trajectory/hops_trajectory.py @@ -1,25 +1,35 @@ +from __future__ import annotations + import copy import os import time as timer import warnings +from pathlib import Path +from typing import Dict, List, Optional, Sequence, Tuple, Union, Any import numpy as np import scipy as sp -import scipy.sparse as sps +import scipy.sparse as sparse +from mesohops.basis.hops_aux import AuxiliaryVector as AuxVec from mesohops.basis.hops_basis import HopsBasis -from mesohops.basis.hops_hierarchy import HopsHierarchy +from mesohops.basis.hops_hierarchy import HopsHierarchy, HIERARCHY_DICT_DEFAULT from mesohops.basis.hops_system import HopsSystem from mesohops.eom.hops_eom import HopsEOM +from mesohops.noise.hops_noise import NOISE_DICT_DEFAULT from mesohops.noise.prepare_functions import prepare_hops_noise from mesohops.storage.hops_storage import HopsStorage from mesohops.util.dynamic_dict import Dict_wDefaults -from mesohops.util.exceptions import LockedException, TrajectoryError, UnsupportedRequest +from mesohops.util.exceptions import ( + LockedException, + TrajectoryError, + UnsupportedRequest, +) from mesohops.util.physical_constants import precision # Constant __title__ = "HOPS" -__author__ = "D. I. G. Bennett, L. Varvelo" -__version__ = "1.2" +__author__ = "D. I. G. B. Raccah, L. Varvelo, J. K. Lynd" +__version__ = "1.6" INTEGRATION_DICT_DEFAULT = { "INTEGRATOR": "RUNGE_KUTTA", @@ -68,17 +78,19 @@ class HopsTrajectory: 'noise2', # Secondary noise trajectory (real) '_z_mem', # Noise memory drift 'noise_param', # Noise parameters (type, seed, etc.) + 'noise2_param', # Noise parameters for purely-real noise 2 (type, seed, etc.) ) def __init__( self, - system_param, - eom_param=None, - noise_param=None, - hierarchy_param=None, - storage_param=None, - integration_param=None, - ): + system_param: dict[str, Any] | str | os.PathLike[str] | Path | None = None, + eom_param: Dict[str, object] | None = None, + noise_param: Dict[str, object] | None = None, + noise2_param: Dict[str, object] | None = None, + hierarchy_param: Dict[str, object] | None = None, + storage_param: Dict[str, object] | None = None, + integration_param: Dict[str, object] | None = None, + ) -> None: """ This class manages four classes: 1. Class: HopsNoise1 (Hierarchy Noise) @@ -165,29 +177,48 @@ def __init__( self._phi = [] self._z_mem = [] self._t = 0 + # Instantiates sub-classes # ----------------------- eom = HopsEOM(eom_param) system = HopsSystem(system_param) hierarchy = HopsHierarchy(hierarchy_param, system.param) + self.noise_param = noise_param - self.basis = HopsBasis(system, hierarchy, eom) - self.storage = HopsStorage(self.basis.eom.param['ADAPTIVE'],storage_param) - self.noise1 = prepare_hops_noise(noise_param, self.basis.system.param) - if "L_NOISE2" in system.param.keys(): - noise_param2 = copy.copy(noise_param) - rand = np.random.RandomState(seed=noise_param['SEED']) - noise_param2['SEED'] = int(rand.randint(0, 2 ** 20, 1)[0]) - self.noise2 = prepare_hops_noise(noise_param2, self.basis.system.param, flag=2) - else: - noise_param2 = { + + # Defaults the second noise to zero + if noise2_param is None: + noise2_param = { "TLEN": noise_param["TLEN"], "TAU": noise_param["TAU"], "MODEL": "ZERO", + "SEED": None, } - if "NOISE_WINDOW" in noise_param.keys(): - noise_param2["NOISE_WINDOW"] = noise_param["NOISE_WINDOW"] - self.noise2 = prepare_hops_noise(noise_param2, self.basis.system.param, flag=1) + # Warn if the user is employing risky practices + else: + if "SEED" in noise_param.keys() and "SEED" in noise2_param.keys(): + if (noise2_param['SEED'] == noise_param['SEED'] and noise_param['SEED'] is + not None): + warnings.warn("Using the same seed for both noise 1 and " + "noise 2 may introduce unphysical correlations between " + "independent noise terms.") + if (noise2_param["TAU"] != noise_param["TAU"] or noise2_param["TLEN"] != + noise_param["TLEN"]): + if "INTERPOLATE" not in noise2_param.keys(): + warnings.warn("Time axes of noise 1 and noise 2 are mismatched. " + "Be cautious when choosing propagation time step to avoid " + "exceptions.") + elif not noise2_param["INTERPOLATE"]: + warnings.warn("Time axes of noise 1 and noise 2 are mismatched. " + "Be cautious when choosing propagation time step to avoid " + "exceptions.") + self.noise2_param = noise2_param + + self.basis = HopsBasis(system, hierarchy, eom) + self.storage = HopsStorage(self.basis.eom.param['ADAPTIVE'],storage_param) + self.noise1 = prepare_hops_noise(self.noise_param, self.basis.system.param) + self.noise2 = prepare_hops_noise(self.noise2_param, self.basis.system.param, + noise_type=2) # Defines integration method # ------------------------- @@ -219,7 +250,11 @@ def __init__( # LOCKING VARIABLE self.__initialized__ = False - def initialize(self, psi_0, timer_checkpoint=None): + def initialize( + self, + psi_0: Sequence[complex] | np.ndarray, + timer_checkpoint: float | None = None, + ) -> None: """ Initializes the trajectory module by ensuring that each sub-component is prepared to begin propagating a trajectory. @@ -270,14 +305,10 @@ def initialize(self, psi_0, timer_checkpoint=None): list_stable_state = self.state_list list_state_new = list( set(self.state_list).union(set(self.static_basis[0]))) - list_add_state = set(list_state_new) - set(list_stable_state) - state_update = (list_state_new, list_stable_state, list_add_state) list_stable_aux = self.auxiliary_list list_aux_new = list( set(self.auxiliary_list).union(set(self.static_basis[1]))) - list_add_aux = set(list_aux_new) - set(list_stable_aux) - aux_update = (list_aux_new, list_stable_aux, list_add_aux) (phi_tmp, dsystem_dt) = self.basis.update_basis( phi_tmp, list_state_new, list_aux_new @@ -304,8 +335,14 @@ def initialize(self, psi_0, timer_checkpoint=None): else: raise LockedException("HopsTrajectory.initialize()") - def make_adaptive(self, delta_a=1e-4, delta_s=1e-4, update_step=1, - f_discard=0.01, adaptive_noise=True): + def make_adaptive( + self, + delta_a: float = 1e-4, + delta_s: float = 1e-4, + update_step: int = 1, + f_discard: float = 0.01, + adaptive_noise: bool = True, + ) -> None: """ Transforms a not-yet-initialized HOPS trajectory from a standard HOPS to an adaptive HOPS approach. @@ -361,7 +398,12 @@ def make_adaptive(self, delta_a=1e-4, delta_s=1e-4, update_step=1, else: raise TrajectoryError("Calling make_adaptive on an initialized trajectory") - def propagate(self, t_advance, tau, timer_checkpoint=None): + def propagate( + self, + t_advance: float, + tau: float, + timer_checkpoint: Optional[float] = None, + ) -> None: """ Performs the integration along fixed time-points. The kind of integration that is performed is controlled by 'step' which was setup in the initialization. @@ -441,7 +483,7 @@ def propagate(self, t_advance, tau, timer_checkpoint=None): step_num = 0 while (set(state_update) != set(self.basis.system.state_list) or set(aux_update) != set(self.basis.hierarchy.auxiliary_list)): - state_update, aux_update, phi = self.inchworm_integrate( + state_update, aux_update, phi, z_mem = self.inchworm_integrate( state_update, aux_update, tau ) if self.early_integrator == 'STATIC_STATE_INCHWORM_HIERARCHY': @@ -481,11 +523,11 @@ def propagate(self, t_advance, tau, timer_checkpoint=None): (phi, self.dsystem_dt) = self.basis.update_basis( phi, state_update, aux_update ) - - self.storage.store_step( - phi_new=phi, aux_list=self.auxiliary_list, state_list=self.state_list, t_new=t, - z_mem_new=self.z_mem - ) + if self.storage.check_storage_time(t): + self.storage.store_step( + phi_new=phi, aux_list=self.auxiliary_list, state_list=self.state_list, t_new=t, + z_mem_new=self.z_mem + ) self.phi = phi self.z_mem = z_mem self.t = t @@ -495,7 +537,7 @@ def propagate(self, t_advance, tau, timer_checkpoint=None): self.storage.metadata["LIST_PROPAGATION_TIME"].append(timer.time() - timer_checkpoint) - def _operator(self, op): + def _operator(self, op: np.ndarray | sparse.spmatrix) -> None: """ Acts an operator on the full hierarchy. Automatically adds all states that become populated to the basis in the adaptive case to avoid uncontrolled error, @@ -534,7 +576,7 @@ def _operator(self, op): self.basis.update_basis(self.phi, self.state_list, self.auxiliary_list) self.reset_early_time_integrator() - def _check_tau_step(self, tau, precision): + def _check_tau_step(self, tau: float, precision: float) -> bool: """ Checks if tau_step is within precision of noise1.param['TAU'] and if tau is greater than or equal to noise1.param['TAU']. @@ -571,7 +613,7 @@ def _check_tau_step(self, tau, precision): return p_bool1 and p_bool2 and s_bool1 and s_bool2 - def normalize(self, phi): + def normalize(self, phi: np.ndarray) -> np.ndarray: """ Re-normalizes the wave function at each step to correct for loss of norm due to finite numerical accuracy. @@ -591,7 +633,12 @@ def normalize(self, phi): else: return phi - def inchworm_integrate(self, list_state_new, list_aux_new, tau): + def inchworm_integrate( + self, + list_state_new: Sequence[int], + list_aux_new: Sequence[AuxVec], + tau: float, + ) -> Tuple[list[int], list[AuxVec], np.ndarray, sparse.sparray]: """ Performs inchworm integration. @@ -617,6 +664,9 @@ def inchworm_integrate(self, list_state_new, list_aux_new, tau): 3. phi : np.array(complex) Full state of the hierarchy normalized if appropriate. + 4. z_mem : np.array(complex) + Noise memory drift terms for the bath [units: cm^-1]. + """ # Incorporate the new states into the state basis state_update = list(set(list_state_new) | set(self.state_list)) @@ -638,9 +688,11 @@ def inchworm_integrate(self, list_state_new, list_aux_new, tau): # Define new basis for step in extended space z_step = self._prepare_zstep(z_mem) (state_update, aux_update) = self.basis.define_basis(phi, tau, z_step) - return state_update, aux_update, phi + return state_update, aux_update, phi, z_mem - def _prepare_zstep(self, z_mem): + def _prepare_zstep( + self, z_mem: sparse.spmatrix | np.ndarray + ) -> List[Union[np.ndarray, sparse.spmatrix]]: """ Constructs a compressed version of the noise terms to be used in the following time step @@ -661,7 +713,9 @@ def _prepare_zstep(self, z_mem): z_rnd2 = self.noise2.get_noise([t], list_absindex_L2)[:, 0] return [z_rnd1, z_rnd2, z_mem] - def construct_noise_correlation_function(self, n_l2, n_traj): + def construct_noise_correlation_function( + self, n_l2: int, n_traj: int + ) -> Tuple[np.ndarray, np.ndarray]: """ Uses correlated noise trajectories to reconstruct the full noise correlation function. @@ -696,15 +750,15 @@ def construct_noise_correlation_function(self, n_l2, n_traj): noise1.get_noise(t_axis, self.basis.system.param["LIST_INDEX_L2_BY_HMODE"])[n_l2, :], mode='full') list_ct1 += result[result.size // 2:] if "L_NOISE2" in self.basis.system.param.keys(): - noise2 = prepare_hops_noise(self.noise_param, self.basis.system.param, - flag=2) + noise2 = prepare_hops_noise(self.noise2_param, + self.basis.system.param, noise_type=2) noise2._prepare_noise() result = np.correlate(noise2.get_noise(t_axis, self.basis.system.param["LIST_INDEX_L2_BY_HMODE"])[n_l2, :], noise2.get_noise(t_axis, self.basis.system.param["LIST_INDEX_L2_BY_HMODE"])[n_l2, :], mode='full') list_ct2 += result[result.size // 2:] return list_ct1 / (n_traj * len(t_axis)), list_ct2 / (n_traj * len(t_axis)) - def reset_early_time_integrator(self): + def reset_early_time_integrator(self) -> None: """ Sets self._early_integrator_time to the current time so that the next use of propagate will make the first self.early_steps early time integrator @@ -716,8 +770,15 @@ def reset_early_time_integrator(self): override=True) self._early_step_counter = 0 - def save_slices(self,dir_path_save, list_key=None, file_header=None, seed=None, - step=1, compress=False): + def save_slices( + self, + dir_path_save: str | os.PathLike, + list_key: Sequence[str] | None = None, + file_header: str | None = None, + seed: int | None = None, + step: int = 1, + compress: bool = False, + ) -> None: """ Saves data to disk. By default it will save all the default keys in hops storage dictionary. It can also save specific keys at user defined time intervals. @@ -789,79 +850,277 @@ def save_slices(self,dir_path_save, list_key=None, file_header=None, seed=None, np.savez(filepath, **dict_sliced_data, allow_pickle=True) print(f"Saved sliced data to {filepath}") + def save_checkpoint(self, filepath: str | os.PathLike, + compress: bool = True, + drop_seed: bool = False) -> None: + """ + Save the current state of the trajectory to a file. This involves saving: + a. the input parameters for HopsSystem, HopsHierarchy, HopsNoise1, HopsNoise2, + HopsEOM, HopsIntegration, and HopsStorage objects. + b. the current full wavefunction (phi), + c. the current state basis (state_list), + d. the current auxiliary basis (aux_list), + e. noise memory (z_mem), + f. the simulation time (t), and + g. the early integration counter + h. all data in hops_storage.data + i. all metadata in hops_storage.metadata + + Parameters + ---------- + 1. filepath : str + Output filename for the checkpoint. + 2. compress : bool, optional + If True the file is gzip-compressed. Default is True. + 3. drop_seed : bool, optional + If True the noise seed will not be saved. Default is False. + + Returns + ------- + None + + """ + # HopsEOM Parameter Dictionary: Only keep the default parameters + eom_param = {k: v for k, v in self.basis.eom.param.items() if k != "ADAPTIVE"} + for key in ["DELTA_A", "DELTA_S", "F_DISCARD"]: + if key in eom_param: + eom_param[key] = float(eom_param[key]) + if eom_param.get("UPDATE_STEP") is None: + eom_param.pop("UPDATE_STEP") + + # HopsSystem Parameter Dictionary: + list_hops_sys_param = ['HAMILTONIAN', 'GW_SYSBATH', 'L_HIER', 'L_NOISE1', 'ALPHA_NOISE1', 'PARAM_NOISE1', + 'L_NOISE2', 'ALPHA_NOISE2', 'PARAM_NOISE2', 'L_LT_CORR', 'PARAM_LT_CORR'] + + # Create a dictionary of HopsTrajectory parameters + params = { + 'system_param': {k: self.basis.system.param[k] for k in list_hops_sys_param if k in self.basis.system.param}, + 'hierarchy_param': {k: self.basis.hierarchy.param[k] for k in HIERARCHY_DICT_DEFAULT if k in self.basis.hierarchy.param}, + 'eom_param': eom_param, + 'noise1_param': {k: self.noise1.param[k] for k in NOISE_DICT_DEFAULT if ((k in self.noise1.param) + and not ((k == 'SEED') + and drop_seed + ) + )}, + 'noise2_param': {k: self.noise2.param[k] for k in NOISE_DICT_DEFAULT if ((k in self.noise2.param) + and not ((k == 'SEED') + and drop_seed + ) + )}, + 'integration_param': self.integration_param, + 'storage_param': self.storage.storage_dic, + } + + # Create a dictionary of indexing/bookkeeping variables + checkpoint = { + 'phi': self.phi, + 'z_mem': self.z_mem, + 't': self.t, + 'state_list': np.array(self.state_list, dtype=int), + 'aux_list': np.array([aux.array_aux_vec for aux in self.auxiliary_list], dtype=object), + 'early_counter': self._early_step_counter, + 'storage_data': self.storage.data, + 'storage_meta': self.storage.metadata, + 'params': params, + } + if (not drop_seed) and (self.noise1.param['SEED'] is None): + warnings.warn( + "Noise1 seed is None but drop_seed is False; " + "this may cause unphysical behavior.", + UserWarning + ) + + if compress: + np.savez_compressed(filepath, **checkpoint, allow_pickle=True) + else: + np.savez(filepath, **checkpoint, allow_pickle=True) + + @classmethod + def load_checkpoint(cls, + filename: str | os.PathLike, + add_seed1: int | str | os.PathLike | np.ndarray | None = None, + add_seed2: int | str | os.PathLike | np.ndarray | None = None, + add_system_param: str | os.PathLike | None = None) -> HopsTrajectory: + """ + Loads a trajectory object from a checkpoint file that has been generated using save_checkpoint(). + + Parameters + ---------- + 1. filename : str or os.PathLike + Path to the ``.npz`` file generated using save_checkpoint(). + + 2. add_seed1 : int, str, os.PathLike, np.ndarray or None + Optional parameter to specify a seed value for Noise1, + this will override the seed value (if any) stored in the + checkpoint file. + + 3. add_seed2 : int, str, os.PathLike, np.ndarray or None + Optional parameter to specify a seed value for Noise2, + this will override the seed value (if any) stored in the + checkpoint file. + + 4. add_system_param : str or os.PathLike or None + Optional parameter to specify a path to a file containing a pre-saved + system parameter file. This will enable directly loading the system + parameter dictionary rather than having to pass it through the constructor. + + Returns + ------- + 1. traj : instance(HopsTrajectory) + Reconstructed trajectory with the state and parameters + restored from the checkpoint. + """ + # Load checkpoint data from the .npz file + data = np.load(filename, allow_pickle=True) + params = data['params'].item() + + # Update noise if needed + if add_seed1 is not None: + params['noise1_param']['SEED'] = add_seed1 + + # Update noise if needed + if add_seed2 is not None: + params['noise2_param']['SEED'] = add_seed2 + + # Update system parameters if needed + if add_system_param is not None: + params['system_param'] = add_system_param + + # Instantiate a new trajectory object with the stored parameters + traj = cls( + params['system_param'], + eom_param=params['eom_param'], + noise_param=params['noise1_param'], + noise2_param=params['noise2_param'], + hierarchy_param=params['hierarchy_param'], + storage_param=params['storage_param'], + integration_param=params['integration_param'], + ) + + # Initialize the trajectory with the stored wave function + psi_0 = np.zeros(traj.basis.system.param['NSTATES'], dtype=np.complex128) + psi_0[data['state_list']] = data['phi'][:data['state_list'].size] + traj.initialize(psi_0) + + # Set the auxiliary list based on the stored data + list_aux = [AuxVec(aux, traj.basis.hierarchy.n_hmodes) for aux in data['aux_list']] + # The next block ensures that if an aux is already in the basis it is re-used, rather than + # being defined again. This is consistent with how aux are defined in the adaptive propagation. + for aux_orig in traj.auxiliary_list: + if aux_orig in list_aux: + index_aux = list_aux.index(aux_orig) + list_aux[index_aux] = aux_orig + + # Update the aux basis + # Because of how the hierarchy class is updated, we need to update the basis one level of the + # hierarchy at a time until the entire hierarchy is updated. (There is a requirement that the + # auxiliaries that are added to the basis are within one step of a previously defined aux.) + for depth in range(1,traj.basis.hierarchy.param["MAXHIER"]+1): + list_aux_depth = [aux for aux in list_aux if aux._sum <= depth] + (traj.phi, traj.dsystem_dt) = traj.basis.update_basis(traj.phi, + data['state_list'], + list_aux_depth) + + # The trajectory has the correct basis. Restore the: state vector, + # noise memory and bookkeeping variables. + traj.phi = data['phi'] + traj.z_mem = data['z_mem'].item() + traj.t = float(data['t']) + traj._early_step_counter = int(data['early_counter']) + + # Restore the storage data from the checkpoint file + traj.storage.data = data['storage_data'].item() + traj.storage.metadata = data['storage_meta'].item() + return traj + + def save_system_parameters(self, filepath: str | os.PathLike) -> None: + """ + Saves the system parameters to a file. + + Parameters + ---------- + 1. filepath : str or os.PathLike + Path to the output file where the system parameters will be saved. + + Returns + ------- + None + """ + self.basis.system.save_dict_param(filepath) + @property - def early_integrator(self): + def early_integrator(self) -> str: return self.integration_param["EARLY_ADAPTIVE_INTEGRATOR"] @property - def integrator(self): + def integrator(self) -> str: return self.integration_param["INTEGRATOR"] @property - def inchworm_cap(self): + def inchworm_cap(self) -> int: return self.integration_param["INCHWORM_CAP"] @property - def effective_noise_integration(self): + def effective_noise_integration(self) -> bool: return self.integration_param["EFFECTIVE_NOISE_INTEGRATION"] @property - def static_basis(self): + def static_basis(self) -> list | np.ndarray | None: return self.integration_param["STATIC_BASIS"] @property - def psi(self): + def psi(self) -> np.ndarray: return self.phi[: self.n_state] @property - def n_hier(self): + def n_hier(self) -> int: return self.basis.n_hier @property - def n_state(self): + def n_state(self) -> int: return self.basis.n_state @property - def state_list(self): + def state_list(self) -> list[int]: return self.basis.system.state_list @property - def auxiliary_list(self): + def auxiliary_list(self) -> list[AuxVec]: return self.basis.hierarchy.auxiliary_list @property - def update_step(self): + def update_step(self) -> int | bool | None: return self.basis.eom.param["UPDATE_STEP"] @property - def phi(self): + def phi(self) -> np.ndarray: return self._phi @phi.setter - def phi(self, phi): + def phi(self, phi: np.ndarray) -> None: self._phi = phi @property - def z_mem(self): + def z_mem(self) -> sparse.spmatrix: return self._z_mem @z_mem.setter - def z_mem(self, z_mem): + def z_mem(self, z_mem: sparse.spmatrix) -> None: self._z_mem = z_mem @property - def t(self): + def t(self) -> float: return self._t @property - def early_steps(self): + def early_steps(self) -> int: return self.integration_param["EARLY_INTEGRATOR_STEPS"] @property - def use_early_integrator(self): + def use_early_integrator(self) -> bool: return self._early_step_counter < self.early_steps @t.setter - def t(self, t): + def t(self, t: float) -> None: self._t = t diff --git a/src/mesohops/util/bath_corr_functions.py b/src/mesohops/util/bath_corr_functions.py index 66d8915..aae2d13 100644 --- a/src/mesohops/util/bath_corr_functions.py +++ b/src/mesohops/util/bath_corr_functions.py @@ -12,41 +12,6 @@ # dictionary. -def bcf_convert_sdl_to_exp(lambda_sdl, gamma_sdl, omega_sdl, temp): - """ - Converts a shifted drude-lorentz spectral density parameters to the exponential - equivalent. - - NOTE: THIS WILL NEED TO BE REPLACED WITH A MORE ROBUST FITTING - ROUTINE SIMILAR TO WHAT EISFELD HAS DONE PREVIOUSLY. - - Parameters - ---------- - 1. lambda_sdl : float - Reorganization energy [units: cm^-1]. - - 2. gamma_sdl : float - Reorganization time scale [units: cm^-1]. - - 3. omega_sdl : float - Vibrational frequency [units: cm^-1]. - - 4. temp : float - Temperature [units: K]. - - Returns - ------- - 1. g_exp : complex - Exponential prefactor [units: cm^-2]. - - 2. w_exp : complex - Exponent [units: cm^-1]. - """ - beta = 1 / (kB * temp) - g_exp = 2 * lambda_sdl / beta - 1j * lambda_sdl * gamma_sdl - w_exp = gamma_sdl - 1j * omega_sdl - - return (g_exp, w_exp) def bcf_convert_dl_ud_to_exp(lambda_dl, gamma_dl, omega_dl, temp): """ @@ -102,7 +67,7 @@ def bcf_convert_dl_ud_to_exp(lambda_dl, gamma_dl, omega_dl, temp): return [g_1, -1j*w_1, g_2, -1j*w_2] -def bcf_convert_dl_to_exp_with_Matsubara(lambda_dl, gamma_dl, temp, k_matsubara): +def bcf_convert_dl_to_exp(lambda_dl, gamma_dl, temp, k_matsubara=0): """ Gives the high temperature mode from the Drude-Lorentz spectral density with a user-selected number of Matsubara frequencies and the corresponding corrections to diff --git a/src/mesohops/util/git_utils.py b/src/mesohops/util/git_utils.py new file mode 100644 index 0000000..7637e02 --- /dev/null +++ b/src/mesohops/util/git_utils.py @@ -0,0 +1,40 @@ +import os +import subprocess + +__title__ = "Git Utilities" +__author__ = "A. Hartzell, Z. W. Freeman" +__version__ = "1.6" + + +def get_git_commit_hash(): + """ + Retrieves the current git commit hash of the repository. + + Returns + ------- + 1. commit_hash : str + The current git commit hash, or the git error message if there was + an error retrieving the hash. + """ + # Get the directory of the current file + current_dir = os.path.dirname(os.path.abspath(__file__)) + + try: + # Run git command to get the current commit hash + result = subprocess.run( + ['git', 'rev-parse', 'HEAD'], + cwd=current_dir, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True, + check=False + ) + + if result.returncode == 0: + # If the command was successful, return the commit hash + return result.stdout.strip() + else: + # If the command failed, return the error message + return result.stderr.strip() + except FileNotFoundError: + return "Git is not installed or not found in PATH." diff --git a/src/mesohops/util/helper_functions.py b/src/mesohops/util/helper_functions.py index 07c37c1..29d6a5e 100644 --- a/src/mesohops/util/helper_functions.py +++ b/src/mesohops/util/helper_functions.py @@ -45,4 +45,22 @@ def index_raveler(ls): 1. list_raveled : list(list(int)) Raveled list of integers [[i], [j], [k]...] """ - return [[i] for i in ls] \ No newline at end of file + return [[i] for i in ls] + +def get_states_from_L2(lop): + """ + Fetches the states that the L operators interacts with. + + Parameters + ---------- + 1. lop : np.array(complex) + L2 operator. + + Returns + ------- + 1. tuple : tuple + Tuple of states that correspond to the specific L operator. + """ + + i_x, i_y = np.nonzero(lop) + return tuple(set(i_x) | set(i_y)) diff --git a/tests/test_adap_basis.py b/tests/test_adap_basis.py index 7a12429..75d5594 100644 --- a/tests/test_adap_basis.py +++ b/tests/test_adap_basis.py @@ -1,6 +1,6 @@ import numpy as np from scipy import sparse -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS from mesohops.util.physical_constants import hbar @@ -33,7 +33,7 @@ def const_gw_sysbath(nsite, e_lambda, gamma, temp, gamma_mark): The diagonal Lindblad operator--a site-projection operator. """ # define exponential parameters - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) # define parameter lists gw_sysbath = [] lop_list = [] diff --git a/tests/test_basis_functions_adaptive.py b/tests/test_basis_functions_adaptive.py index 97854f4..d8e49f3 100644 --- a/tests/test_basis_functions_adaptive.py +++ b/tests/test_basis_functions_adaptive.py @@ -134,7 +134,7 @@ def get_error_term_hier_flux_up(list_k_vec, list_w_bymode, list_lop_bymode, P2_p # E[m,k] = |w_m * (k_m+1)|^2 * \sum_d |(L_m @ \psi_k)[d]|^2 E_mode_aux = np.abs(L2_lop) ** 2 @ np.abs(P2_phi[:, aux_index]) ** 2 E_mode_aux *= np.abs((1 + k_mode) * w_mode) ** 2 - E2_err[mode_index, aux_index] = np.sum(E_mode_aux) + E2_err[mode_index, aux_index] = float(np.sum(E_mode_aux)) return E2_err / hbar2 @@ -192,12 +192,17 @@ def get_error_term_state_flux_up(list_k_vec, list_w_bymode, list_lop_bymode, P2_ # s] * \psi_k[s]|^2 for d in range(n_dest): dest = dest_state_list[d] - err_element = np.abs((1 + k_mode) * w_mode) ** 2 * np.abs( - L2_lop[dest, state] * P2_phi[state_index, - aux_index]) ** 2 + err_element = float( + np.abs((1 + k_mode) * w_mode) ** 2 + * np.abs( + L2_lop[dest, state] + * P2_phi[state_index, aux_index] + ) + ** 2 + ) E_state_aux += err_element list_known_err_by_d[d] += err_element - E2_err[state_index, aux_index] = E_state_aux + E2_err[state_index, aux_index] = float(E_state_aux) return E2_err / hbar2, list_known_err_by_d / hbar2 @@ -559,12 +564,15 @@ def get_error_term_hier_flux_down(list_gw_bymode, list_lop_bymode, list_l_exp_by for aux_index in range(n_aux): E_mode_aux = np.abs(L2_lop)**2 @ np.abs(P2_phi[:, aux_index]) ** 2 E_mode_aux *= np.abs(gw_mode)**2 - E2_err[mode_index, aux_index] = np.sum(E_mode_aux) + E2_err[mode_index, aux_index] = float(np.sum(E_mode_aux)) if len(list_s_not_in_d) > 0: for s in list_rel_s_not_in_d: - E2_err[mode_index, aux_index] += np.abs(gw_mode)**2 * np.abs( - list_l_exp_bymode[mode_index]) **2 * np.abs(P2_phi[s, - aux_index]) ** 2 + val = ( + np.abs(gw_mode) ** 2 + * np.abs(list_l_exp_bymode[mode_index]) ** 2 + * np.abs(P2_phi[s, aux_index]) ** 2 + ) + E2_err[mode_index, aux_index] += float(np.sum(val)) return E2_err/hbar2 def get_error_term_state_flux_down(list_gw_bymode, list_lop_bymode, list_l_exp_bymode, @@ -629,17 +637,23 @@ def get_error_term_state_flux_down(list_gw_bymode, list_lop_bymode, list_l_exp_b list_l_exp_bymode[mode_index]) for d in range(n_dest): dest = dest_state_list[d] - err_element = (np.abs(list_gw_bymode[mode_index] * - L2_lop[dest, state_list[state_index]] * - P2_phi[state_index, aux_index])**2) + err_element = np.abs( + list_gw_bymode[mode_index] + * L2_lop[dest, state_list[state_index]] + * P2_phi[state_index, aux_index] + ) ** 2 + err_element = float(np.sum(err_element)) E_state_aux += err_element list_known_err_by_d[d] += err_element if len(list_s_not_in_d) > 0: if state_index in list_rel_s_not_in_d: - E_state_aux += (np.abs(list_gw_bymode[mode_index] * - list_l_exp_bymode[mode_index] * - P2_phi[state_index, aux_index])**2) - E2_err[state_index, aux_index] = E_state_aux + err_add = np.abs( + list_gw_bymode[mode_index] + * list_l_exp_bymode[mode_index] + * P2_phi[state_index, aux_index] + ) ** 2 + E_state_aux += float(np.sum(err_add)) + E2_err[state_index, aux_index] = float(E_state_aux) return E2_err/hbar2, list_known_err_by_d/hbar2 def get_X2_exp_mode_state(list_lop, list_states, list_exp_lop): diff --git a/tests/test_bath_corr_functions.py b/tests/test_bath_corr_functions.py index ffbbce6..ff2399b 100644 --- a/tests/test_bath_corr_functions.py +++ b/tests/test_bath_corr_functions.py @@ -1,25 +1,9 @@ from mesohops.util.bath_corr_functions import * from mesohops.util.physical_constants import kB -def test_bcf_convert_to_sdl(): +def test_bcf_convert_dl_to_exp(): """ - Tests that bcf_convert_to_sdl returns the correct exponential in the correct format. - """ - e_lambda = 100 - gamma = 20 - omega = 33 - temp = 300 - bcf_overdamped = bcf_convert_sdl_to_exp(e_lambda, gamma, 0, temp) - overdamped_analytic = (complex(2*e_lambda*temp*kB - 1j*e_lambda*gamma), - complex(gamma)) - assert overdamped_analytic == bcf_overdamped - bcf_underdamped = bcf_convert_sdl_to_exp(e_lambda, gamma, omega, temp) - underdamped_analytic = (overdamped_analytic[0], gamma-1j*omega) - assert underdamped_analytic == bcf_underdamped - -def test_bcf_convert_dl_to_exp_with_Matsubara(): - """ - Tests that bcf_convert_dl_to_exp_with_Matsubara returns the correct list of + Tests that bcf_convert_dl_to_exp returns the correct list of exponentials in the correct format, and that it approaches the analytic limit of the hyperbolic cotangent real portion of the low-temperature mode. """ @@ -29,10 +13,12 @@ def test_bcf_convert_dl_to_exp_with_Matsubara(): mats_const = np.pi*2*kB*temp # Tests that the 0-Matsubara mode limit reproduces the pure high temperature # approximation - assert bcf_convert_dl_to_exp_with_Matsubara(e_lambda, gamma, temp, 0) == \ - list(bcf_convert_sdl_to_exp(e_lambda, gamma, 0, temp)) + overdamped_analytic = [complex(2 * e_lambda * temp * kB - 1j * e_lambda * gamma), + complex(gamma)] + assert (bcf_convert_dl_to_exp(e_lambda, gamma, temp, 0) == + overdamped_analytic) # Tests that the function returns a list with 2 entries per mode - kmats_10000 = bcf_convert_dl_to_exp_with_Matsubara(e_lambda, gamma, temp, 10000) + kmats_10000 = bcf_convert_dl_to_exp(e_lambda, gamma, temp, 10000) assert len(kmats_10000) == 20002 # Test that Matsubara modes are correct for n in range(3,len(kmats_10000),2): diff --git a/tests/test_checkpoint.py b/tests/test_checkpoint.py new file mode 100644 index 0000000..68d324a --- /dev/null +++ b/tests/test_checkpoint.py @@ -0,0 +1,1197 @@ +import numpy as np +import pytest +import warnings +from scipy import sparse +from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS +from mesohops.trajectory.exp_noise import bcf_exp +from mesohops.storage.hops_storage import HopsStorage +from mesohops.basis.hops_aux import AuxiliaryVector +from .utils import compare_dictionaries + +# List of hierarchy properties that depend on the previous time step +list_hierarchy_properties_path_dependent = [ + '_new_aux_index_conn_by_mode', # New auxiliary index connections by mode + '_new_aux_id_conn_by_mode', # New auxiliary ID connections by mode + '_stable_aux_id_conn_by_mode', # Stable auxiliary ID connections by mode + '_auxiliary_list', # List of current auxiliary vectors (main storage) + '__previous_auxiliary_list', # Auxiliary list from previous step + '__list_aux_stable', # Auxiliaries stable between steps + '__list_aux_remove', # Auxiliaries to remove in update + '__list_aux_add', # Auxiliaries to add in update + '__previous_list_auxstable_index', # Indices of stable auxiliaries from previous step + '_previous_list_modes_in_use', # Modes in use from previous step +] + +list_hierarchy_properties_obj = ['system'] + +list_system_properties_path_dependent = [ + '__previous_state_list', # Previous state list (for adaptive updates) + '__list_add_state', # States to add in update + '__list_stable_state', # States stable between updates + '_list_boundary_state', # States coupled to basis by Hamiltonian + '__list_absindex_new_state_modes', # New state mode indices (absolute) +] +list_system_properties_obj = [] + +list_mode_properties_path_dependent= [] +list_mode_properties_obj = ['system', 'hierarchy'] + +list_aux_properties_path_dependent = [] +list_aux_properties_obj = [] + + +def get_private(obj, name: str): + if name.startswith('__'): + cls = obj.__class__ + mangled = f"_{cls.__name__}{name}" + return getattr(obj, mangled) + else: + return getattr(obj, name) + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def make_hops_nonadaptive(): + """Factory fixture for creating and initialising a non-adaptive HOPS object.""" + # Parameters mirrored from test_hops_trajectory + noise_param = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + } + + noise2_param = { + "SEED": 1010101, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + } + + loperator = np.zeros([2, 2, 2], dtype=np.float64) + loperator[0, 0, 0] = 1.0 + loperator[1, 1, 1] = 1.0 + + + sys_param = { + "HAMILTONIAN": np.array([[0, 10.0], [10.0, 0]], dtype=np.float64), + "GW_SYSBATH": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]], + "L_HIER": [loperator[0], loperator[0], loperator[1], loperator[1]], + "L_NOISE1": [loperator[0], loperator[0], loperator[1], loperator[1]], + "ALPHA_NOISE1": bcf_exp, + "PARAM_NOISE1": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]], + } + + sys_param_with_noise2 = sys_param.copy() + sys_param_with_noise2.update({ + "L_NOISE2": [loperator[0], loperator[0], loperator[1], loperator[1]], + "ALPHA_NOISE2": bcf_exp, + "PARAM_NOISE2": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]], + }) + + system_ltc_params = { + "PARAM_LT_CORR": [250.0 / 1000.0, 250.0 / 2000.0], + "L_LT_CORR": [loperator[0], loperator[1]] + } + + hier_param = {"MAXHIER": 3} + + eom_param = {"TIME_DEPENDENCE": False, "EQUATION_OF_MOTION": "NORMALIZED NONLINEAR"} + + integrator_param = { + "INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 5, + 'INCHWORM_CAP': 5, + 'STATIC_BASIS': None, + 'EFFECTIVE_NOISE_INTEGRATION': False, + } + + psi_0 = [1.0 + 0.0j, 0.0 + 0.0j] + + def _factory(sys_p=[sys_param, sys_param_with_noise2], + system_ltc_params = system_ltc_params, + psi0=psi_0, + noise_p=noise_param, + hier_p=hier_param, + eom_p=eom_param, + integ_p=integrator_param, + storage_param=None, + flag_noise2=False, + noise2_p=noise2_param, + flag_lt_corr=True): + if flag_noise2: + sys_p = sys_p[1] + else: + sys_p = sys_p[0] + noise2_p = None + + if flag_lt_corr: + sys_p.update(system_ltc_params) + + hops = HOPS( + sys_p, + noise_param=noise_p, + noise2_param=noise2_p, + hierarchy_param=hier_p, + eom_param=eom_p, + integration_param=integ_p, + storage_param=storage_param, + ) + hops.initialize(psi0) + return hops + + return _factory + +@pytest.fixture +def make_hops_adaptive(): + nsite = 10 + coupling = 50 + loperator = np.zeros([nsite, nsite, nsite], dtype=np.float64) + for i in range(nsite): + loperator[i, i, i] = 1.0 + + gw_sysbath = [] + lop_list = [] + for i in range(nsite): + gw_sysbath.extend([[20.0, 20.0], [10.0, 10.0]]) + lop_list.extend([loperator[i], loperator[i]]) + + hamiltonian = np.zeros([nsite, nsite], dtype=np.float64) + for i in range(nsite - 1): + hamiltonian[i, i + 1] = coupling + hamiltonian[i + 1, i] = coupling + + sys_param_adapt = { + "HAMILTONIAN": hamiltonian, + "GW_SYSBATH": gw_sysbath, + "L_HIER": lop_list, + "L_NOISE1": lop_list, + "ALPHA_NOISE1": bcf_exp, + "PARAM_NOISE1": gw_sysbath, + } + + psi0 = np.zeros(nsite, dtype=np.complex128) + psi0[0] = 1.0 + 0.0j + + noise_param = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + } + + hier_param = {"MAXHIER": 3} + + eom_param = { + "TIME_DEPENDENCE": False, + "EQUATION_OF_MOTION": "NORMALIZED NONLINEAR", + "DELTA_A": 1e-3, + "DELTA_S": 1e-3, + "UPDATE_STEP": 10, + "ADAPTIVE_S": True, + "ADAPTIVE_H": True, + } + + integrator_param = { + "INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 5, + 'INCHWORM_CAP': 5, + 'STATIC_BASIS': None, + 'EFFECTIVE_NOISE_INTEGRATION': False, + } + + def _factory(sys_p=sys_param_adapt, + psi0=psi0, + noise_p=noise_param, + hier_p=hier_param, + eom_p=eom_param, + integ_p=integrator_param, + storage_param=None, + ): + hops = HOPS( + sys_p, + noise_param=noise_p, + hierarchy_param=hier_p, + eom_param=eom_p, + integration_param=integ_p, + storage_param=storage_param, + ) + hops.initialize(psi0) + return hops + + return _factory + +def test_compare_dictionaries(tmp_path, make_hops_nonadaptive): + """ + This test checks that compare_dictionaries() correctly manages the comparison + between dictionaries with a variety of different data types and nesting. + + It tests: + 1. Dictionaries with different keys + 2. Dictionaries with the same keys but different values + 3. Dictionaries with the same keys but different value types + 4. Dictionaries with nested structures that differ + 5. Dictionaries with arrays or lists that differ + 6. Edge cases like empty dictionaries, None values, etc. + + Specific test cases include: + - Test 1: Identical dictionaries should compare equal + Verifies that comparing a dictionary with itself works correctly. + + - Test 2: Dictionaries with different keys should raise an error + Compares non-adaptive and adaptive storage dictionaries which have different keys. + + - Test 3: Dictionaries with the same keys but different values + Modifies a time value in a copy of a dictionary and verifies that comparison fails. + + - Test 4: Dictionaries with the same keys but different value types + Tests that a list and a numpy array with the same values are considered equal. + + - Test 5: Dictionaries with arrays that differ slightly + Verifies that small differences (within default tolerance) are ignored. + + - Test 6: Dictionaries with arrays that differ significantly + Confirms that larger differences cause the comparison to fail. + + - Test 7: Dictionaries with nested structures + Tests comparison of dictionaries containing nested dictionaries with differences. + + - Test 8: Empty dictionaries + Verifies that empty dictionaries compare equal. + + - Test 9: None values + Tests comparison of dictionaries containing None values. + + - Test 10: Sparse matrices + Verifies comparison of dictionaries containing scipy sparse matrices. + + - Test 11: Lists of different lengths + Tests that lists with different lengths are detected as different. + + - Test 12: Object arrays + Tests comparison of numpy object arrays containing dictionaries. + + - Additional tests with actual HOPS objects: + Verifies that checkpoint data matches the original and that modifications + to the loaded data are detected correctly. + """ + # Create a HopsStorage object with non-adaptive mode + storage_non_adaptive = HopsStorage(False, { + 'psi_traj': True, + 't_axis': True, + 'z_mem': True, + }) + + # Create a HopsStorage object with adaptive mode + storage_adaptive = HopsStorage(True, { + 'psi_traj': True, + 't_axis': True, + 'z_mem': True, + 'aux_list': True, + 'state_list': True, + 'list_nhier': True, + 'list_nstate': True, + 'list_aux_norm': True, + }) + + # Set dimensions for both storage objects + storage_non_adaptive.n_dim = 5 + storage_adaptive.n_dim = 5 + + # Create test data + phi_new = np.array([1.0, 2.0, 3.0, 4.0, 5.0], dtype=np.complex128) + aux_list = [AuxiliaryVector([], 4), AuxiliaryVector([[0, 1]], 4)] + state_list = [0, 1, 2, 3, 4] + t_new = 1.0 + z_mem_new = np.array([1.0, 2.0, 3.0], dtype=np.complex128) + + # Store data in both storage objects + storage_non_adaptive.store_step( + phi_new=phi_new, + aux_list=aux_list, + state_list=state_list, + t_new=t_new, + z_mem_new=z_mem_new + ) + + storage_adaptive.store_step( + phi_new=phi_new, + aux_list=aux_list, + state_list=state_list, + t_new=t_new, + z_mem_new=z_mem_new + ) + + # Create reference dictionaries for testing + dict1 = storage_non_adaptive.data.copy() + dict2 = storage_adaptive.data.copy() + + # Test 1: Dictionaries with the same content should compare equal + compare_dictionaries(dict1, dict1) + compare_dictionaries(dict2, dict2) + print("Test 1 passed: Identical dictionaries compare equal") + + # Test 2: Dictionaries with different keys should raise an error + with pytest.raises(AssertionError, match="Dictionaries must have the same keys"): + compare_dictionaries(dict1, dict2) + print("Test 2 passed: Dictionaries with different keys raise the expected error") + + # Test 3: Dictionaries with the same keys but different values + dict3 = dict1.copy() + dict3['t_axis'] = [2.0] # Change the time value + with pytest.raises(AssertionError): + compare_dictionaries(dict1, dict3) + print("Test 3 passed: Dictionaries with different values raise an error") + + # Test 4: Dictionaries with the same keys but different value types + dict4 = dict1.copy() + dict4['t_axis'] = np.array([1.0]) # Change list to numpy array + compare_dictionaries(dict1, dict4) + print("Test 4 passed: Dictionaries with compatible types compare equal") + + # Test 5: Dictionaries with arrays that differ slightly + dict5 = dict1.copy() + dict5['psi_traj'] = [np.array(dict1['psi_traj'][0]) + 1e-10] # Small difference + compare_dictionaries(dict1, dict5) + print("Test 5 passed: Arrays with small differences compare equal within tolerance") + + # Test 6: Dictionaries with arrays that differ significantly + dict6 = dict1.copy() + dict6['psi_traj'] = [np.array(dict1['psi_traj'][0]) + 0.1] # Larger difference + with pytest.raises(AssertionError): + compare_dictionaries(dict1, dict6) + print("Test 6 passed: Arrays with significant differences raise an error") + + # Test 7: Dictionaries with nested structures + dict7 = {'nested': {'a': 1, 'b': 2}, 'array': np.array([1, 2, 3])} + dict8 = {'nested': {'a': 1, 'b': 3}, 'array': np.array([1, 2, 3])} + with pytest.raises(AssertionError): + compare_dictionaries(dict7, dict8) + print("Test 7 passed: Nested dictionaries with differences raise an error") + + # Test 8: Edge case - empty dictionaries + compare_dictionaries({}, {}) + print("Test 8 passed: Empty dictionaries compare equal") + + # Test 9: Edge case - None values + dict9 = {'a': None} + dict10 = {'a': None} + dict11 = {'a': 0} + compare_dictionaries(dict9, dict10) + with pytest.raises(AssertionError): + compare_dictionaries(dict9, dict11) + print("Test 9 passed: Dictionaries with None values handled correctly") + + # Test 10: Sparse matrices + sparse_mat1 = sparse.csr_array(([1, 2, 3], ([0, 1, 2], [0, 1, 2])), shape=(3, 3)) + sparse_mat2 = sparse.csr_array(([1, 2, 3], ([0, 1, 2], [0, 1, 2])), shape=(3, 3)) + sparse_mat3 = sparse.csr_array(([1, 2, 4], ([0, 1, 2], [0, 1, 2])), shape=(3, 3)) + + dict12 = {'sparse': sparse_mat1} + dict13 = {'sparse': sparse_mat2} + dict14 = {'sparse': sparse_mat3} + + compare_dictionaries(dict12, dict13) + with pytest.raises(AssertionError): + compare_dictionaries(dict12, dict14) + print("Test 10 passed: Dictionaries with sparse matrices handled correctly") + + # Test 11: Lists of different lengths + dict15 = {'list': [1, 2, 3]} + dict16 = {'list': [1, 2]} + with pytest.raises(AssertionError): + compare_dictionaries(dict15, dict16) + print("Test 11 passed: Lists of different lengths raise an error") + + + # Test 12: Object arrays + obj_array1 = np.array([{'a': 1}, {'b': 2}], dtype=object) + obj_array2 = np.array([{'a': 1}, {'b': 2}], dtype=object) + obj_array3 = np.array([{'a': 1}, {'b': 3}], dtype=object) + + dict17 = {'obj_array': obj_array1} + dict18 = {'obj_array': obj_array2} + dict19 = {'obj_array': obj_array3} + + compare_dictionaries(dict17, dict18) + with pytest.raises(AssertionError): + compare_dictionaries(dict17, dict19) + print("Test 12 passed: Object arrays handled correctly") + + + # Test with actual HOPS objects + hops = make_hops_nonadaptive() + hops.propagate(100.0, 1.0) + + ckpt_path = tmp_path / "traj.npz" + hops.save_checkpoint(str(ckpt_path)) + + storage_mid = {k: list(v) if isinstance(v, list) else v for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + for key in storage_mid: + np.testing.assert_array_equal(hops_loaded.storage.data[key], storage_mid[key]) + print("Test with HOPS objects passed: Checkpoint data matches original") + + + # Modify the loaded data to test failure cases + if 't_axis' in hops_loaded.storage.data: + original_t_axis = hops_loaded.storage.data['t_axis'].copy() + hops_loaded.storage.data['t_axis'] = [t + 0.1 for t in hops_loaded.storage.data['t_axis']] + + with pytest.raises(AssertionError): + for key in storage_mid: + np.testing.assert_array_equal(hops_loaded.storage.data[key], storage_mid[key]) + print("Test with modified HOPS data passed: Modified data detected correctly") + + # Restore the original data + hops_loaded.storage.data['t_axis'] = original_t_axis + + + +def test_checkpoint_nonadaptive(tmp_path, make_hops_nonadaptive): + """Checkpoints for a standard, non-adaptive trajectory.""" + + hops = make_hops_nonadaptive() + hops.propagate(100.0, 1.0) + + ckpt_path = tmp_path / "traj.npz" + hops.save_checkpoint(str(ckpt_path)) + + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + hops.propagate(100.0, 1.0) + phi_final = hops.phi.copy() + storage_final = hops.storage.data + t_final = hops.t + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + + # Test storage_mid values + for key in storage_mid: + np.testing.assert_array_equal(hops_loaded.storage.data[key], storage_mid[key]) + + hops_loaded.propagate(100.0, 1.0) + + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + for key in storage_final: + np.testing.assert_array_equal(hops_loaded.storage.data[key], storage_final[key]) + + +def test_checkpoint_early_time_integration(tmp_path, make_hops_nonadaptive): + """Ensures checkpoints work during the early-time integration phase.""" + + # Create a trajectory and propagate for fewer steps than EARLY_INTEGRATOR_STEPS + hops = make_hops_nonadaptive(integ_p={"EARLY_ADAPTIVE_INTEGRATOR": "INCH_WORM", + "EARLY_INTEGRATOR_STEPS": 50}) + hops.propagate(50.0, 2.0) + + assert hops.use_early_integrator # Early integrator should still be active + early_counter_mid = hops._early_step_counter + + ckpt_path = tmp_path / "traj_early.npz" + hops.save_checkpoint(str(ckpt_path)) + + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Continue propagation beyond the early integration window + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + t_final = hops.t + early_counter_final = hops._early_step_counter + storage_final = hops.storage.data + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Confirm the loaded trajectory matches the checkpoint state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + assert hops_loaded._early_step_counter == early_counter_mid + assert hops_loaded.use_early_integrator + compare_dictionaries(storage_mid, hops_loaded.storage.data) + + # Continue propagation from the loaded state + hops_loaded.propagate(100., 2.0) + + # Final states should match the original trajectory + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + assert hops_loaded._early_step_counter == early_counter_final + compare_dictionaries(storage_final, hops_loaded.storage.data) + + +def test_checkpoint_without_lt_corr(tmp_path, make_hops_nonadaptive): + """Tests checkpointing without low temperature corrections.""" + # Initialize HOPS with low temperature corrections + hops = make_hops_nonadaptive(flag_lt_corr=False) + + # Propagate the system + hops.propagate(100.0, 1.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_lt_corr.npz" + hops.save_checkpoint(str(ckpt_path)) + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + lt_corr_param_mid = hops.basis.system.list_lt_corr_param.copy() + + # Continue propagation + hops.propagate(100.0, 1.0) + phi_final = hops.phi.copy() + t_final = hops.t + storage_final = hops.storage.data + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + + # Test storage_mid values + for key in storage_mid: + np.testing.assert_array_equal(hops_loaded.storage.data[key], storage_mid[key]) + + # Continue propagation from the loaded state + hops_loaded.propagate(100.0, 1.0) + + # Verify the final state matches the expected final state + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + for key in storage_final: + np.testing.assert_array_equal(hops_loaded.storage.data[key], storage_final[key]) + +def test_checkpoint_uncompressed(tmp_path, make_hops_nonadaptive): + """Verifies checkpoints saved without compression.""" + + hops = make_hops_nonadaptive() + hops.propagate(50.0, 1.0) + + ckpt_path = tmp_path / "traj_uncompressed.npz" + hops.save_checkpoint(str(ckpt_path), compress=False) + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + np.testing.assert_allclose(hops_loaded.phi, hops.phi) + assert hops_loaded.t == hops.t + +def test_checkpoint_missing_directory(tmp_path, make_hops_nonadaptive): + """Ensures saving to a missing directory raises an error.""" + + hops = make_hops_nonadaptive() + ckpt_path = tmp_path / "missing" / "traj.npz" + with pytest.raises(FileNotFoundError): + hops.save_checkpoint(str(ckpt_path)) + +def test_checkpoint_adaptive_storage(tmp_path, make_hops_adaptive): + """Check storage recovery with adaptive propagation.""" + hops = make_hops_adaptive( + storage_param={ + 'phi_traj': True, + 'psi_traj': True, + 't_axis': True, + 'z_mem': True, + 'aux_list': True, + 'state_list': True, + 'list_nhier': True, + 'list_nstate': True, + 'list_aux_norm': True, + }, + ) + + # Propagate with adaptive calculation + hops.propagate(50.0, 2.0) + + ckpt_path = tmp_path / "traj_adaptive.npz" + hops.save_checkpoint(str(ckpt_path)) + + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + storage_final = hops.storage.data + t_final = hops.t + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + + compare_dictionaries(storage_mid, hops_loaded.storage.data) + assert len(hops_loaded.auxiliary_list[-1]) == len(hops.auxiliary_list[-1]) + + hops_loaded.propagate(100.0, 2.0) + + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + compare_dictionaries(storage_final, hops_loaded.storage.data) + + +def test_checkpoint_adaptive_hierarchy(tmp_path, make_hops_adaptive): + """Ensures hierarchy objects are restored correctly in adaptive runs.""" + + hops = make_hops_adaptive() + + # Propagate with adaptive calculation + hops.propagate(50.0, 2.0) + list_hierarchy_properties = [str for str in hops.basis.hierarchy.__class__.__slots__] + list_hierarchy_properties_mid = [str for str in list_hierarchy_properties + if not (str in list_hierarchy_properties_path_dependent) and + not (str in list_hierarchy_properties_obj)] + list_hierarchy_properties_fin = [str for str in list_hierarchy_properties + if not (str in list_hierarchy_properties_obj)] + original_hierarchy_mid = {prop: get_private(hops.basis.hierarchy, prop) + for prop in list_hierarchy_properties_mid} + + ckpt_path = tmp_path / "traj_adaptive.npz" + hops.save_checkpoint(str(ckpt_path)) + hops.propagate(100.0, 2.0) + original_hierarchy_fin = {prop: get_private(hops.basis.hierarchy, prop) + for prop in list_hierarchy_properties_fin} + + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + loaded_hierarchy_mid = {prop: get_private(hops_loaded.basis.hierarchy, prop) + for prop in list_hierarchy_properties_mid} + hops_loaded.propagate(100.0, 2.0) + loaded_hierarchy_fin = {prop: get_private(hops_loaded.basis.hierarchy, prop) + for prop in list_hierarchy_properties_fin} + + # Compare the hierarchy properties + compare_dictionaries(original_hierarchy_mid, loaded_hierarchy_mid) + compare_dictionaries(original_hierarchy_fin, loaded_hierarchy_fin) + + +def test_checkpoint_adaptive_modes(tmp_path, make_hops_adaptive): + """Check that mode information is preserved across checkpoints.""" + hops = make_hops_adaptive() + + # Propagate with adaptive calculation + hops.propagate(50.0, 2.0) + + list_mode_properties = [str for str in hops.basis.mode.__class__.__slots__] + list_mode_properties_mid = [str for str in list_mode_properties + if not (str in list_mode_properties_path_dependent) and + not (str in list_mode_properties_obj)] + list_mode_properties_fin = [str for str in list_mode_properties + if not (str in list_mode_properties_obj)] + + original_mode_mid = {prop: get_private(hops.basis.mode, prop) + for prop in list_mode_properties_mid} + + ckpt_path = tmp_path / "traj_adaptive.npz" + hops.save_checkpoint(str(ckpt_path)) + hops.propagate(100.0, 2.0) + original_mode_fin = {prop: get_private(hops.basis.mode, prop) + for prop in list_mode_properties_fin} + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + loaded_mode_mid = {prop: get_private(hops_loaded.basis.mode, prop) + for prop in list_mode_properties_mid} + hops_loaded.propagate(100.0, 2.0) + loaded_mode_fin = {prop: get_private(hops_loaded.basis.mode, prop) + for prop in list_mode_properties_fin} + + # Compare the hierarchy properties + compare_dictionaries(original_mode_mid, loaded_mode_mid) + compare_dictionaries(original_mode_fin, loaded_mode_fin) + + +def test_checkpoint_adaptive_system(tmp_path, make_hops_adaptive): + """Verifies system data is reconstructed with adaptive parameters.""" + hops = make_hops_adaptive() + + ckpt_path = tmp_path / "traj_adaptive.npz" + + # Propagate with adaptive calculation + hops.propagate(50.0, 2.0) + + list_system_properties = [str for str in hops.basis.system.__class__.__slots__] + list_system_properties_mid = [str for str in list_system_properties + if not (str in list_system_properties_path_dependent) and + not (str in list_system_properties_obj)] + list_system_properties_fin = [str for str in list_system_properties + if not (str in list_system_properties_obj)] + + hops.save_checkpoint(str(ckpt_path)) + orig_system_mid = {prop: get_private(hops.basis.system, prop) + for prop in list_system_properties_mid} + hops.propagate(100.0, 2.0) + orig_system_fin = {prop: get_private(hops.basis.system, prop) + for prop in list_system_properties_fin} + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + load_system_mid = {prop: get_private(hops_loaded.basis.system, prop) + for prop in list_system_properties_mid} + hops_loaded.propagate(100.0, 2.0) + load_system_fin = {prop: get_private(hops_loaded.basis.system, prop) + for prop in list_system_properties_fin} + + compare_dictionaries(orig_system_mid, load_system_mid) + compare_dictionaries(orig_system_fin, load_system_fin) + +def test_checkpoint_adaptive_listaux(tmp_path, make_hops_adaptive): + """Verifies auxiliary lists survive checkpointing in adaptive mode.""" + hops = make_hops_adaptive() + + ckpt_path = tmp_path / "traj_adaptive.npz" + + # Propagate with adaptive calculation + hops.propagate(50.0, 2.0) + list_aux_properties = [str for str in hops.auxiliary_list[0].__class__.__slots__] + list_aux_properties_mid = [str for str in list_aux_properties + if not (str in list_aux_properties_path_dependent) and + not (str in list_aux_properties_obj)] + list_aux_properties_fin = [str for str in list_aux_properties + if not (str in list_aux_properties_obj)] + hops.save_checkpoint(str(ckpt_path)) + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + for aux_orig, aux_load in zip(hops.auxiliary_list, hops_loaded.auxiliary_list): + dict_aux_orig = {key: get_private(aux_orig, key) for key in list_aux_properties_mid} + dict_aux_load = {key: get_private(aux_load, key) for key in list_aux_properties_mid} + compare_dictionaries(dict_aux_orig, dict_aux_load) + + hops.propagate(100.0, 2.0) + hops_loaded.propagate(100.0, 2.0) + for aux_orig, aux_load in zip(hops.auxiliary_list, hops_loaded.auxiliary_list): + dict_aux_orig = {key: get_private(aux_orig, key) for key in list_aux_properties_fin} + dict_aux_load = {key: get_private(aux_load, key) for key in list_aux_properties_fin} + compare_dictionaries(dict_aux_orig, dict_aux_load) + + +def test_checkpoint_orphan_aux(tmp_path, make_hops_adaptive): + """Auxiliary without parents should load and propagate correctly.""" + # Create a trajectory and propagate briefly to generate a checkpoint + hops = make_hops_adaptive() + hops.propagate(20.0, 2.0) + + + list_aux_properties = [str for str in hops.auxiliary_list[0].__class__.__slots__] + list_aux_properties_mid = [str for str in list_aux_properties + if not (str in list_aux_properties_path_dependent) and + not (str in list_aux_properties_obj)] + list_aux_properties_fin = [str for str in list_aux_properties + if not (str in list_aux_properties_obj)] + + + list_hierarchy_properties = [str for str in hops.basis.hierarchy.__class__.__slots__] + list_hierarchy_properties_mid = [str for str in list_hierarchy_properties + if not (str in list_hierarchy_properties_path_dependent) and + not (str in list_hierarchy_properties_obj)] + list_hierarchy_properties_fin = [str for str in list_hierarchy_properties + if not (str in list_hierarchy_properties_obj)] + + # Construct an auxiliary list composed entirely of orphans + list_aux = [aux for aux in hops.auxiliary_list if (aux._sum == 2 or aux._sum == 0)] + phi_tmp, dsystem_dt = hops.basis.update_basis(hops.phi, hops.state_list, list_aux) + hops.phi = phi_tmp + hops.dsystem_dt = dsystem_dt + + ckpt_path = tmp_path / "orig.npz" + hops.save_checkpoint(str(ckpt_path)) + original_hierarchy_mid = {prop: get_private(hops.basis.hierarchy, prop) + for prop in list_hierarchy_properties_mid} + + # Loading the modified checkpoint should succeed even though the + # auxiliary list lacks connections + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + loaded_hierarchy_mid = {prop: get_private(hops_loaded.basis.hierarchy, prop) + for prop in list_hierarchy_properties_mid} + for aux_orig, aux_load in zip(hops.auxiliary_list, hops_loaded.auxiliary_list): + dict_aux_orig = {key: get_private(aux_orig, key) for key in list_aux_properties_mid} + dict_aux_load = {key: get_private(aux_load, key) for key in list_aux_properties_mid} + compare_dictionaries(dict_aux_orig, dict_aux_load) + + compare_dictionaries(hops.storage.data, hops_loaded.storage.data) + compare_dictionaries(original_hierarchy_mid, loaded_hierarchy_mid) + + + # Propagate further to ensure the trajectory is functional + hops_loaded.propagate(50.0, 2.0) + hops.propagate(50.0, 2.0) + loaded_hierarchy_fin = {prop: get_private(hops_loaded.basis.hierarchy, prop) + for prop in list_hierarchy_properties_fin} + original_hierarchy_fin = {prop: get_private(hops.basis.hierarchy, prop) + for prop in list_hierarchy_properties_fin} + + compare_dictionaries(hops.storage.data, hops_loaded.storage.data) + compare_dictionaries(original_hierarchy_fin, loaded_hierarchy_fin) + for aux_orig, aux_load in zip(hops.auxiliary_list, hops_loaded.auxiliary_list): + dict_aux_orig = {key: get_private(aux_orig, key) for key in list_aux_properties_fin} + dict_aux_load = {key: get_private(aux_load, key) for key in list_aux_properties_fin} + compare_dictionaries(dict_aux_orig, dict_aux_load) + + + +def test_checkpoint_with_noise2(tmp_path, make_hops_nonadaptive): + # Initialize HOPS with updated system parameters + hops = make_hops_nonadaptive(flag_noise2=True) + + assert hops.noise2.param["SEED"] == 1010101 + + # Propagate the system + hops.propagate(100.0, 2.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_noise2.npz" + hops.save_checkpoint(str(ckpt_path)) + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + + # Continue propagation + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + storage_final = hops.storage.data + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + assert hops_loaded.noise2.param["SEED"] == 1010101 + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + compare_dictionaries(storage_mid, hops_loaded.storage.data) + + # Continue propagation from the loaded state + hops_loaded.propagate(100.0, 2.0) + + # Verify the final state matches the expected final state + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + compare_dictionaries(storage_final, hops_loaded.storage.data) + + +def test_checkpoint_interpolate(tmp_path, make_hops_nonadaptive): + """Tests checkpointing with INTERPOLATE flag set to True.""" + # Create custom noise parameters with INTERPOLATE=True + noise_param = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + "INTERPOLATE": True, + } + + # Initialize HOPS with INTERPOLATE=True + hops = make_hops_nonadaptive(noise_p=noise_param) + + # Propagate the system + hops.propagate(100.0, 2.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_interpolate.npz" + hops.save_checkpoint(str(ckpt_path)) + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Continue propagation + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + t_final = hops.t + storage_final = hops.storage.data + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + compare_dictionaries(storage_mid, hops_loaded.storage.data) + + # Verify the INTERPOLATE flag was correctly restored + assert hops_loaded.noise1.param["INTERPOLATE"] == True + + # Continue propagation from the loaded state + hops_loaded.propagate(100.0, 2.0) + + # Verify the final state matches the expected final state + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + compare_dictionaries(storage_final, hops_loaded.storage.data) + +def test_checkpoint_adaptive_noise(tmp_path, make_hops_nonadaptive): + """Tests checkpointing with ADAPTIVE noise flag set to True.""" + # Create custom noise parameters with ADAPTIVE=True + noise_param = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + "ADAPTIVE": True, + } + + # Initialize HOPS with ADAPTIVE=True + hops = make_hops_nonadaptive(noise_p=noise_param) + + # Propagate the system + hops.propagate(100.0, 2.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_adaptive_noise.npz" + hops.save_checkpoint(str(ckpt_path)) + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Continue propagation + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + t_final = hops.t + storage_final = hops.storage.data + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + compare_dictionaries(storage_mid, hops_loaded.storage.data) + + # Verify the ADAPTIVE flag was correctly restored + assert hops_loaded.noise1.param["ADAPTIVE"] == True + + # Continue propagation from the loaded state + hops_loaded.propagate(100.0, 2.0) + + # Verify the final state matches the expected final state + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + compare_dictionaries(storage_final, hops_loaded.storage.data) + +def test_checkpoint_store_raw_noise(tmp_path, make_hops_nonadaptive): + """Tests checkpointing with STORE_RAW_NOISE flag set to True.""" + # Create custom noise parameters with STORE_RAW_NOISE=True + noise_param = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + "STORE_RAW_NOISE": True, + } + + # Initialize HOPS with STORE_RAW_NOISE=True + hops = make_hops_nonadaptive(noise_p=noise_param) + + # Propagate the system + hops.propagate(100.0, 2.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_store_raw_noise.npz" + hops.save_checkpoint(str(ckpt_path)) + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Continue propagation + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + t_final = hops.t + storage_final = hops.storage.data + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + compare_dictionaries(storage_mid, hops_loaded.storage.data) + + # Verify the STORE_RAW_NOISE flag was correctly restored + assert hops_loaded.noise1.param["STORE_RAW_NOISE"] == True + + # Continue propagation from the loaded state + hops_loaded.propagate(100.0, 2.0) + + # Verify the final state matches the expected final state + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + compare_dictionaries(storage_final, hops_loaded.storage.data) + +def test_checkpoint_noise_window(tmp_path, make_hops_nonadaptive): + """Tests checkpointing with NOISE_WINDOW parameter set.""" + # Create custom noise parameters with NOISE_WINDOW=100.0 + noise_param = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + "NOISE_WINDOW": 100.0, + } + + # Initialize HOPS with NOISE_WINDOW=100.0 + hops = make_hops_nonadaptive(noise_p=noise_param) + + # Propagate the system + hops.propagate(100.0, 2.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_noise_window.npz" + hops.save_checkpoint(str(ckpt_path)) + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Continue propagation + hops.propagate(100.0, 2.0) + phi_final = hops.phi.copy() + t_final = hops.t + storage_final = hops.storage.data + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + compare_dictionaries(storage_mid, hops_loaded.storage.data) + + # Verify the NOISE_WINDOW parameter was correctly restored + assert hops_loaded.noise1.param["NOISE_WINDOW"] == 100.0 + + # Continue propagation from the loaded state + hops_loaded.propagate(100.0, 2.0) + + # Verify the final state matches the expected final state + np.testing.assert_allclose(hops_loaded.phi, phi_final, atol=1e-100) + assert hops_loaded.t == t_final + compare_dictionaries(storage_final, hops_loaded.storage.data) + +def test_checkpoint_seed_none(tmp_path, make_hops_nonadaptive): + """Tests checkpointing with SEED=None.""" + # Create custom noise parameters with SEED=None + noise_param = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + } + + # Initialize HOPS with SEED=None + hops = make_hops_nonadaptive(noise_p=noise_param) + + # Propagate the system + hops.propagate(100.0, 2.0) + + # Save the checkpoint + ckpt_path = tmp_path / "traj_seed_none.npz" + hops.save_checkpoint(str(ckpt_path)) + + # Save the mid-propagation state + phi_mid = hops.phi.copy() + t_mid = hops.t + storage_mid = {k: list(v) if isinstance(v, list) else v + for k, v in hops.storage.data.items() + if k != 'ADAPTIVE'} + + # Load the checkpoint + hops_loaded = HOPS.load_checkpoint(str(ckpt_path)) + + # Verify the loaded state matches the mid-propagation state + np.testing.assert_allclose(hops_loaded.phi, phi_mid) + assert hops_loaded.t == t_mid + compare_dictionaries(storage_mid, hops_loaded.storage.data) + assert hops_loaded.noise1.param['SEED'] is None + +def test_checkpoint_drop_seed_array(tmp_path, make_hops_nonadaptive): + hops = make_hops_nonadaptive() + ckpt_keep = tmp_path / "keep.npz" + ckpt_drop = tmp_path / "drop.npz" + hops.save_checkpoint(str(ckpt_keep), drop_seed=False) + hops.save_checkpoint(str(ckpt_drop), drop_seed=True) + + params_keep = np.load(ckpt_keep, allow_pickle=True)["params"].item() + params_drop = np.load(ckpt_drop, allow_pickle=True)["params"].item() + + assert np.array_equal(params_keep["noise1_param"]["SEED"], hops.noise1.param["SEED"]) + assert np.array_equal(params_keep["noise2_param"]["SEED"], hops.noise2.param["SEED"]) + assert "SEED" not in params_drop["noise1_param"] + assert "SEED" not in params_drop["noise2_param"] + + +def test_load_checkpoint_add_seed(tmp_path, make_hops_nonadaptive, capsys): + hops = make_hops_nonadaptive(flag_noise2=True) + assert hops.noise2.param["SEED"] == 1010101 + ckpt_path = tmp_path / "add_seed.npz" + hops.save_checkpoint(str(ckpt_path), drop_seed=True) + + new_seed1 = 7 + new_seed2 = 1010108 + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path), add_seed1=new_seed1) + assert hops_loaded.noise1.param["SEED"] == new_seed1 + assert hops_loaded.noise2.param["SEED"] == None + + hops_loaded = HOPS.load_checkpoint(str(ckpt_path), add_seed1=new_seed1, + add_seed2=new_seed2) + assert hops_loaded.noise1.param["SEED"] == 7 + assert hops_loaded.noise2.param["SEED"] == 1010108 + + with warnings.catch_warnings(record=True) as w: + hops_loaded = HOPS.load_checkpoint(str(ckpt_path), add_seed1=new_seed1, + add_seed2=new_seed1) + assert any("Using the same seed for both noise 1 and" in str( + warning.message) for warning in w) + + +def test_noise1_seed_warning(tmp_path, make_hops_nonadaptive): + """Test that a UserWarning is raised when noise1 seed is None and drop_seed is False.""" + # Create custom noise parameters with SEED=None + noise_param = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 500.0, + "TAU": 0.5, + } + + # Initialize HOPS with SEED=None + hops = make_hops_nonadaptive(noise_p=noise_param) + + # Propagate the system + hops.propagate(50.0, 2.0) + + ckpt_path = tmp_path / "checkpoint.npz" + with pytest.warns(UserWarning): + hops.save_checkpoint(ckpt_path, drop_seed=False) diff --git a/tests/test_dimer_of_dimers.py b/tests/test_dimer_of_dimers.py index 63989bd..4593fe8 100644 --- a/tests/test_dimer_of_dimers.py +++ b/tests/test_dimer_of_dimers.py @@ -1,11 +1,12 @@ import os + import numpy as np import scipy as sp -from scipy import sparse +from mesohops.eom.eom_hops_ksuper import _permute_aux_by_matrix from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS -from mesohops.eom.eom_hops_ksuper import _permute_aux_by_matrix -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp +from scipy import sparse __title__ = "Test of eom_integrator_rk_nonlin_norm" __author__ = "D. I. G. Bennett" @@ -28,7 +29,7 @@ e_lambda = 20.0 gamma = 50.0 temp = 140.0 -(g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) +(g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([4, 4, 4], dtype=np.float64) gw_sysbath = [] @@ -206,13 +207,7 @@ def test_integration_variables(): flag_pass = True for key in var_list_desk.keys(): - try: - np.testing.assert_allclose(var_list_desk[key], var_list_lap[key]) - print(key) - except: - print(key, " <-- ERROR!") - flag_pass = False - assert flag_pass == True + np.testing.assert_allclose(var_list_desk[key], var_list_lap[key], err_msg=f"{key} <-- ERROR!") def test_alpha(): diff --git a/tests/test_dyadic_spectra.py b/tests/test_dyadic_spectra.py index 73610d5..99461a9 100644 --- a/tests/test_dyadic_spectra.py +++ b/tests/test_dyadic_spectra.py @@ -1,19 +1,19 @@ -import pytest import numpy as np +import pytest +from scipy import sparse + from mesohops.trajectory.dyadic_spectra import DyadicSpectra as DHOPS from mesohops.trajectory.dyadic_spectra import (prepare_spectroscopy_input_dict, - prepare_chromophore_input_dict, - prepare_convergence_parameter_dict) + prepare_chromophore_input_dict, + prepare_convergence_parameter_dict) from mesohops.util.bath_corr_functions import ishizaki_decomposition_bcf_dl -from scipy import sparse - def test_DyadicSpectra(): """ Tests the DyadicSpectra class for properly unpacking input dictionaries, and ensures the Hamiltonian is the proper shape. """ - ## Spectroscopy input dictionary + # Spectroscopy input dictionary seed = 10 spectrum_type = "FLUORESCENCE" propagation_time_dict = {"t_2": 2.0, "t_3": 3.0} @@ -24,7 +24,7 @@ def test_DyadicSpectra(): propagation_time_dict, field_dict, site_dict) - ## Chromophore input dictionary + # Chromophore input dictionary M2_mu_ge = np.array([np.array([0.5, 0.2, 0.1]), np.array([0.5, 0.2, 0.1])]) H2_sys_hamiltonian = np.zeros((3, 3), dtype=np.complex128) H2_sys_hamiltonian[1:, 1:] = np.array([[0, -100], [-100, 0]]) @@ -46,13 +46,13 @@ def test_DyadicSpectra(): bath_dict2) # Case 3: list_modes + static_filter_list (Markovian) - static_filter_list = ['Markovian', [True, False]] + static_filter_list = [['Markovian', [True, False]]] bath_dict3 = {"list_lop": list_lop, "list_modes": list_modes, "static_filter_list": static_filter_list} chromophore_dict_3 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict3) - ## Convergence parameter dictionaries + # Convergence parameter dictionaries convergence_dict_float_dt = prepare_convergence_parameter_dict(t_step=0.1, max_hier=12, delta_a=1e-10, @@ -66,7 +66,7 @@ def test_DyadicSpectra(): set_update_step=1, set_f_discard=0.5) - ## Test input dictionary unpacking + # Test input dictionary unpacking # Case 1 + float dt dhops_1a = DHOPS(spectroscopy_dict, chromophore_dict_1, convergence_dict_float_dt, @@ -142,27 +142,27 @@ def test_DyadicSpectra(): assert np.allclose(dhops_3.ltc_param, chromophore_dict_3["ltc_param"]) for a, b in zip(dhops_3.lop_list_ltc, chromophore_dict_3["lop_list_ltc"]): assert np.allclose(a.toarray(), b.toarray()) - assert dhops_3.static_filter_list == ['Markovian', [True, False]] + assert dhops_3.static_filter_list == [['Markovian', [True, False, True, False]]] - ## Test Hamiltonian shape compatibility with number of states + # Test Hamiltonian shape compatibility with number of states H2_sys_hamiltonian_wrongshape = np.zeros((4, 4), dtype=np.complex128) bath_dict_wrongshape = {"list_lop": list_lop, "list_modes": list_modes} - chromophore_dict_wrongshape = prepare_chromophore_input_dict(M2_mu_ge, - H2_sys_hamiltonian_wrongshape, - bath_dict_wrongshape) - - try: + chromophore_dict_wrongshape = ( + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian_wrongshape, + bath_dict_wrongshape)) + + with pytest.raises(ValueError, + match='H2_sys_hamiltonian must be \\(\\(n_chrom \\+ 1\\) x ' + '\\(n_chrom \\+ 1\\)\\) to account for each chromophore and' + ' the ground state.'): DHOPS(spectroscopy_dict, chromophore_dict_wrongshape, convergence_dict_float_dt, seed) - except ValueError as excinfo: - if ('H2_sys_hamiltonian must be ((n_chrom + 1) x (n_chrom + 1)) to account' - ' for each chromophore and the ground state.') not in str(excinfo): - pytest.fail() # Test "INITIALIZATION_TIME" is greater than 0 - dhops_4 = DHOPS(spectroscopy_dict, chromophore_dict_1, convergence_dict_float_dt, seed) + dhops_4 = DHOPS(spectroscopy_dict, chromophore_dict_1, convergence_dict_float_dt, + seed) dhops_4.calculate_spectrum() assert dhops_4.storage.metadata["INITIALIZATION_TIME"] > 0 @@ -171,11 +171,13 @@ def test_DyadicSpectra(): assert len(dhops_4.storage.metadata["LIST_PROPAGATION_TIME"]) == 2 # Test "LIST_PROPAGATION_TIME" is correct length (Absorption) - spectroscopy_dict_abs = prepare_spectroscopy_input_dict("ABSORPTION", - {"t_1": 1.0}, - {"E_1": np.array([0, 0, 1])}, - {"list_ket_sites": np.array([1, 2])}) - dhops_5 = DHOPS(spectroscopy_dict_abs, chromophore_dict_1, convergence_dict_float_dt, seed) + spectroscopy_dict_abs = ( + prepare_spectroscopy_input_dict("ABSORPTION", + {"t_1": 1.0}, + {"E_1": np.array([0, 0, 1])}, + {"list_ket_sites": np.array([1, 2])})) + dhops_5 = DHOPS(spectroscopy_dict_abs, chromophore_dict_1, + convergence_dict_float_dt, seed) dhops_5.calculate_spectrum() assert len(dhops_5.storage.metadata["LIST_PROPAGATION_TIME"]) == 1 @@ -186,7 +188,7 @@ def test_initialize(capsys): Tests the initialization of the DyadicTrajectory class, and ensures that the class is properly initialized in both non-adaptive and adaptive cases. """ - ## Spectroscopy input dictionary + # Spectroscopy input dictionary seed = 10 spectrum_type = "ABSORPTION" propagation_time_dict = {"t_1": 1.0} @@ -197,7 +199,7 @@ def test_initialize(capsys): propagation_time_dict, field_dict, site_dict) - ## Chromophore input dictionary + # Chromophore input dictionary M2_mu_ge = np.array([np.array([0.5, 0.2, 0.1]), np.array([0.5, 0.2, 0.1])]) H2_sys_hamiltonian = np.zeros((3, 3), dtype=np.complex128) H2_sys_hamiltonian[1:, 1:] = np.array([[0, -100], [-100, 0]]) @@ -210,23 +212,23 @@ def test_initialize(capsys): chromophore_dict = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict) - ## Convergence parameter dictionaries + # Convergence parameter dictionaries convergence_dict = prepare_convergence_parameter_dict(t_step=0.1, max_hier=12, delta_a=0, delta_s=0) - ## Testing initialized property works upon initialization + # Testing initialized property works upon initialization dhops = DHOPS(spectroscopy_dict, chromophore_dict, convergence_dict, seed) assert dhops.__initialized__ is False dhops.initialize() assert dhops.__initialized__ is True - ## Testing that multiple calls to initialize() triggers a warning + # Testing that multiple calls to initialize() triggers a warning dhops.initialize() out, err = capsys.readouterr() assert "WARNING: DyadicTrajectory has already been initialized." in out.strip() - ## Testing delta_a/delta_s greater than 0 causes the trajectory to be ran adaptively + # Testing delta_a/delta_s greater than 0 causes the trajectory to be ran adaptively convergence_dict_adaptive = prepare_convergence_parameter_dict(t_step=0.1, max_hier=12, delta_a=1e-10, @@ -234,7 +236,8 @@ def test_initialize(capsys): set_update_step=2, set_f_discard=0.5) - dhops_adaptive = DHOPS(spectroscopy_dict, chromophore_dict, convergence_dict_adaptive, seed) + dhops_adaptive = DHOPS(spectroscopy_dict, chromophore_dict, + convergence_dict_adaptive, seed) dhops_adaptive.initialize() # Adaptive @@ -253,7 +256,7 @@ def test_hilb_operator(): Tests the _hilb_operator method of the DyadicTrajectory class, ensuring that the Hilbert raising and lowering operators are properly constructed. """ - ## Spectroscopy input dictionary + # Spectroscopy input dictionary seed = 10 spectrum_type = "ABSORPTION" propagation_time_dict = {"t_1": 1.0} @@ -264,7 +267,7 @@ def test_hilb_operator(): propagation_time_dict, field_dict, site_dict) - ## Chromophore input dictionary + # Chromophore input dictionary M2_mu_ge = np.array([np.array([0.5, 0.2, 0.1]), np.array([0.5, 0.2, 0.1])]) H2_sys_hamiltonian = np.zeros((3, 3), dtype=np.complex128) H2_sys_hamiltonian[1:, 1:] = np.array([[0, -100], [-100, 0]]) @@ -277,13 +280,13 @@ def test_hilb_operator(): chromophore_dict = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict) - ## Convergence parameter dictionaries + # Convergence parameter dictionaries convergence_dict = prepare_convergence_parameter_dict(t_step=0.1, max_hier=12, delta_a=0, delta_s=0) dhops = DHOPS(spectroscopy_dict, chromophore_dict, convergence_dict, seed) - ## Testing that the Hilbert raising and lowering operators are properly constructed + # Testing that the Hilbert raising and lowering operators are properly constructed # Note that the value of 1.7 comes from the dot product of the field and the dipole dense_raise = np.zeros((3, 3), dtype=np.float64) @@ -293,24 +296,24 @@ def test_hilb_operator(): dense_lower[0, np.array([1, 2])] = 1.7 assert np.allclose(dhops._hilb_operator("raise", np.array([2, 3, 1]), - dhops.list_ket_sites).toarray(), dense_raise) + dhops.list_ket_sites).toarray(), + dense_raise) assert np.allclose(dhops._hilb_operator("lower", np.array([2, 3, 1]), - dhops.list_ket_sites).toarray(), dense_lower) + dhops.list_ket_sites).toarray(), + dense_lower) - ## Testing that the method raises an error if not given a valid action_type - try: + # Testing that the method raises an error if not given a valid action_type + with pytest.raises(ValueError, match="action_type must be either 'raise' or " + "'lower'."): dhops._hilb_operator("cha_cha_slide", np.array([2, 3, 1]), dhops.list_ket_sites) - except ValueError as excinfo: - if "action_type must be either 'raise' or 'lower'." not in str(excinfo): - pytest.fail() def test_final_dyad_operator(): """ Tests the _final_dyad_operator method of the DyadicTrajectory class, ensuring that the final dyad operator is properly constructed, and the time index is properly set. """ - ## Spectroscopy input dictionary + # Spectroscopy input dictionary seed = 10 spectrum_type = "FLUORESCENCE" propagation_time_dict = {"t_2": 2.0, "t_3": 3.0} @@ -321,7 +324,7 @@ def test_final_dyad_operator(): propagation_time_dict, field_dict, site_dict) - ## Chromophore input dictionary + # Chromophore input dictionary M2_mu_ge = np.array([np.array([0.5, 0.2, 0.1]), np.array([0.5, 0.2, 0.1])]) H2_sys_hamiltonian = np.zeros((3, 3), dtype=np.complex128) H2_sys_hamiltonian[1:, 1:] = np.array([[0, -100], [-100, 0]]) @@ -334,13 +337,13 @@ def test_final_dyad_operator(): chromophore_dict = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict) - ## Convergence parameter dictionaries + # Convergence parameter dictionaries convergence_dict = prepare_convergence_parameter_dict(t_step=0.1, max_hier=12, delta_a=1e-10, delta_s=1e-10, set_update_step=1, set_f_discard=0.5) - ## Testing that the final dyad operator is properly constructed + # Testing that the final dyad operator is properly constructed dhops = DHOPS(spectroscopy_dict, chromophore_dict, convergence_dict, seed) dyadic_f_op = np.zeros((6, 6), dtype=np.float64) @@ -349,7 +352,7 @@ def test_final_dyad_operator(): assert np.allclose(dhops._final_dyad_operator()[0].toarray(), dyadic_f_op) - ## Testing that the time index is properly set + # Testing that the time index is properly set assert dhops._final_dyad_operator()[1] == 20 def test_prepare_spectroscopy_input_dict(capsys): @@ -412,14 +415,11 @@ def test_prepare_spectroscopy_input_dict(capsys): # Testing site definition errors # Case 1: list_ket_sites not defined - try: + with pytest.raises(ValueError, match='list_ket_sites must be defined.'): prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, propagation_time_dict={"t_1": t1}, field_dict={"E_1": E1}, site_dict={}) - except ValueError as excinfo: - if 'list_ket_sites must be defined.' not in str(excinfo): - pytest.fail() # Case 2: list_ket_sites not a numpy array ket_list = prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, @@ -436,61 +436,49 @@ def test_prepare_spectroscopy_input_dict(capsys): assert np.allclose(ket_list["list_ket_sites"], ket_array["list_ket_sites"]) # Case 3: sites indexed from 0 - try: + with pytest.raises(ValueError, match='Ket and Bra sites must be indexed starting ' + 'from 1.'): prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, propagation_time_dict={"t_1": t1}, field_dict={"E_1": E1}, site_dict={ "list_ket_sites": ket_sites_index_issue}) - except ValueError as excinfo: - if 'Ket and Bra sites must be indexed starting from 1.' not in str(excinfo): - pytest.fail() # Testing field input formatting # Case 1: Field not a numpy array - try: + with pytest.raises(ValueError, match='All field entries should be numpy arrays.'): prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, propagation_time_dict={"t_1": t1}, field_dict={"E_1": E1_list}, site_dict={"list_ket_sites": ket_sites}) - except ValueError as excinfo: - if 'All field entries should be numpy arrays.' not in str(excinfo): - pytest.fail() # Case 2: Field not a numpy array with exactly 3 entries - try: + with pytest.raises(ValueError, + match='All field entries should be numpy arrays with exactly 3 ' + 'entries.'): prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, propagation_time_dict={"t_1": t1}, field_dict={"E_1": E1_wrong_length}, site_dict={"list_ket_sites": ket_sites}) - except ValueError as excinfo: - if ('All field entries should be numpy arrays with exactly 3 entries.' not in - str(excinfo)): - pytest.fail() # Testing under-defined absorption input # Case 1: t_1 not defined - try: + with pytest.raises(ValueError, + match='Propagation time after first field interaction \\(t_1\\) ' + 'must be defined as > 0 for absorption.'): prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, propagation_time_dict={}, field_dict={"E_1": E1}, site_dict={"list_ket_sites": ket_sites}) - except ValueError as excinfo: - if ('Propagation time after first field interaction (t_1) must be defined as > ' - '0 for absorption.') not in str(excinfo): - pytest.fail() # Case 2: E_1 not defined - try: + with pytest.raises(ValueError, match='E_1 must be defined for absorption.'): prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, propagation_time_dict={"t_1": t1}, field_dict={}, site_dict={"list_ket_sites": ket_sites}) - except ValueError as excinfo: - if 'E_1 must be defined for absorption.' not in str(excinfo): - pytest.fail() # Testing over-defined absorption input @@ -501,7 +489,7 @@ def test_prepare_spectroscopy_input_dict(capsys): site_dict={"list_ket_sites": ket_sites}) out, err = capsys.readouterr() assert out.strip() == ('WARNING: Only t_1 is necessary for absorption. ' - 'Setting all other propagation times to zero.') + 'Setting all other propagation times to zero.') # Case 2: field_dict contains too many inputs prepare_spectroscopy_input_dict(spectrum_type=absorption_spectrum_type, @@ -510,19 +498,17 @@ def test_prepare_spectroscopy_input_dict(capsys): site_dict={"list_ket_sites": ket_sites}) out, err = capsys.readouterr() assert out.strip() == ('WARNING: Only E_1 is necessary for absorption. E_sig is ' - 'set to E_1. All other field definitions will be discarded') + 'set to E_1. All other field definitions will be ' + 'discarded') # Testing under-defined fluorescence input # Case 1: list_bra_sites not defined - try: + with pytest.raises(ValueError, match='list_bra_sites must be defined for fluorescence.'): prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, propagation_time_dict={"t_2": t2, "t_3": t3}, field_dict={"E_1": E1, "E_sig": Esig}, site_dict={"list_ket_sites": ket_sites}) - except ValueError as excinfo: - if 'list_bra_sites must be defined for fluorescence.' not in str(excinfo): - pytest.fail() # Case 2: list_bra_sites not a numpy array bra_list = prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, @@ -545,38 +531,33 @@ def test_prepare_spectroscopy_input_dict(capsys): # already tested prior to determining the spectrum type. # Case 3: t_2 or t_3 not defined - try: + with pytest.raises(ValueError, + match='Propagation times after second and third field ' + 'interactions \\(t_2, t_3\\) must be defined as > 0 for ' + 'fluorescence.'): prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, propagation_time_dict={"t_3": t3}, field_dict={"E_1": E1, "E_sig": Esig}, site_dict={"list_ket_sites": ket_sites, "list_bra_sites": bra_sites}) - except ValueError as excinfo: - if ('Propagation times after second and third field interactions (t_2, t_3) ' - 'must be defined as > 0 for fluorescence.') not in str(excinfo): - pytest.fail() - try: + with pytest.raises(ValueError, + match='Propagation times after second and third field ' + 'interactions \\(t_2, t_3\\) must be defined as > 0 for ' + 'fluorescence.'): prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, propagation_time_dict={"t_2": t2}, field_dict={"E_1": E1, "E_sig": Esig}, site_dict={"list_ket_sites": ket_sites, "list_bra_sites": bra_sites}) - except ValueError as excinfo: - if ('Propagation times after second and third field interactions (t_2, t_3) ' - 'must be defined as > 0 for fluorescence.') not in str(excinfo): - pytest.fail() # Case 4: E_1 not defined - try: + with pytest.raises(ValueError, match='E_1 must be defined for fluorescence.'): prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, propagation_time_dict={"t_2": t2, "t_3": t3}, field_dict={"E_sig": Esig}, site_dict={"list_ket_sites": ket_sites, "list_bra_sites": bra_sites}) - except ValueError as excinfo: - if 'E_1 must be defined for fluorescence.' not in str(excinfo): - pytest.fail() # Case 5: E_sig not defined prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, @@ -586,7 +567,7 @@ def test_prepare_spectroscopy_input_dict(capsys): "list_bra_sites": bra_sites}) out, err = capsys.readouterr() assert out.strip() == ('WARNING: E_sig is not defined. Setting E_sig to default, ' - '[0, 0, 1].') + '[0, 0, 1].') # Testing over-defined fluorescence input @@ -599,7 +580,7 @@ def test_prepare_spectroscopy_input_dict(capsys): "list_bra_sites": bra_sites}) out, err = capsys.readouterr() assert out.strip() == ('WARNING: Only t_2 and t_3 are necessary for fluorescence. ' - 'Setting all other propagation times to zero.') + 'Setting all other propagation times to zero.') # Case 2: field_dict contains too many inputs prepare_spectroscopy_input_dict(spectrum_type=fluorescence_spectrum_type, @@ -609,18 +590,15 @@ def test_prepare_spectroscopy_input_dict(capsys): "list_bra_sites": bra_sites}) out, err = capsys.readouterr() assert out.strip() == ('WARNING: Only E_1 and E_sig are necessary for fluorescence.' - ' All other field definitions will be discarded.') + ' All other field definitions will be discarded.') # Testing incorrect spectrum_type input - try: + with pytest.raises(ValueError, match='spectrum_type must be one of the following:'): prepare_spectroscopy_input_dict(spectrum_type=bad_spectrum_type, propagation_time_dict={"t_2": t2, "t_3": t3}, field_dict={"E_1": E1, "E_sig": Esig}, site_dict={"list_ket_sites": ket_sites, "list_bra_sites": bra_sites}) - except ValueError as excinfo: - if 'spectrum_type must be one of the following:' not in str(excinfo): - pytest.fail() def test_prepare_chromophore_input_dict(): @@ -630,42 +608,53 @@ def test_prepare_chromophore_input_dict(): """ kmax2 = 4 + kmax2_negative = -2 M2_mu_ge = np.array([[0.5, 0.2, 0.1], [0.5, 0.2, 0.1]]) H2_sys_hamiltonian = np.zeros((3, 3), dtype=np.complex128) H2_sys_hamiltonian[1:, 1:] = np.array([[0, -100], [-100, 0]]) list_lop = [sparse.coo_matrix(([1], ([1], [2])), shape=(3, 3)), sparse.coo_matrix(([1], ([2], [1])), shape=(3, 3))] - lop_list_noise_nmodes_LTC = [sparse.coo_matrix(([1], ([1], [2])), shape=(3, 3)), + lop_list_noise = [sparse.coo_matrix(([1], ([1], [2])), shape=(3, 3)), sparse.coo_matrix(([1], ([1], [2])), shape=(3, 3)), sparse.coo_matrix(([1], ([2], [1])), shape=(3, 3)), sparse.coo_matrix(([1], ([2], [1])), shape=(3, 3))] list_modes = ishizaki_decomposition_bcf_dl(35, 50, 295, 0) - list_modes_by_bath = [ishizaki_decomposition_bcf_dl(40, 60, 290, 0), ishizaki_decomposition_bcf_dl(35, 50, 295, 0)] + list_modes_by_bath = [ishizaki_decomposition_bcf_dl(40, 60, 290, 0), + ishizaki_decomposition_bcf_dl(35, 50, 295, 0)] nmodes_LTC = 1 - static_filter_list_markovian = ['Markovian', [True, False]] - static_filter_list_triangular = ['Triangular', [[True, False], kmax2]] - static_filter_list_longedge = ['LongEdge', [[True, False], kmax2]] - - ## Testing proper output - - # Case 1: list_lop, list_modes, nmodes_LTC, and static_filter_list - bath_dict_1 = {"list_lop": list_lop, "list_modes": list_modes, "nmodes_LTC": nmodes_LTC, "static_filter_list": static_filter_list_longedge} - chromophore_dict_1 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict_1) + static_filter_list_lm = [['Markovian', [True, False]], + ['Triangular', [[True, False], kmax2]], + ['LongEdge', [[True, False], kmax2]]] + static_filter_list_lmbb = [['Markovian', [True, False, True, False]], + ['Triangular', [[True, False, True, False], kmax2]], + ['LongEdge', [[True, False, True, False], kmax2]]] + + # Testing proper output + + # Case 1: list_lop, list_modes, and static_filter_list_lm + bath_dict_1 = {"list_lop": list_lop, "list_modes": list_modes, + "static_filter_list": static_filter_list_lm} + chromophore_dict_1 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + bath_dict_1) assert np.allclose(chromophore_dict_1["M2_mu_ge"], M2_mu_ge) assert chromophore_dict_1["n_chromophore"] == 2 assert np.allclose(chromophore_dict_1["H2_sys_hamiltonian"], H2_sys_hamiltonian) - # Note: when nmodes_LTC is set to 1, lop_list_hier=list_lop=lop_list_ltc for a dimer - for a, b in zip(chromophore_dict_1["lop_list_hier"], list_lop): + # Note: when nmodes_LTC is not set, lop_list_hier=lop_list_noise and lop_list_ltc=[] + for a, b in zip(chromophore_dict_1["lop_list_hier"], lop_list_noise): assert np.allclose(a.toarray(), b.toarray()) - for a, b in zip(chromophore_dict_1["lop_list_ltc"], list_lop): + for a, b in zip(chromophore_dict_1["lop_list_ltc"], []): assert np.allclose(a.toarray(), b.toarray()) - for a, b in zip(chromophore_dict_1["lop_list_noise"], lop_list_noise_nmodes_LTC): + for a, b in zip(chromophore_dict_1["lop_list_noise"], lop_list_noise): assert np.allclose(a.toarray(), b.toarray()) assert np.allclose(chromophore_dict_1["gw_sysbath_hier"], [(list_modes[0], - list_modes[1]), - (list_modes[0], - list_modes[1])]) + list_modes[1]), + (list_modes[2], + list_modes[3]), + (list_modes[0], + list_modes[1]), + (list_modes[2], + list_modes[3])]) assert np.allclose(chromophore_dict_1["gw_sysbath_noise"], [(list_modes[0], list_modes[1]), (list_modes[2], @@ -674,66 +663,57 @@ def test_prepare_chromophore_input_dict(): list_modes[1]), (list_modes[2], list_modes[3])]) - assert np.allclose(chromophore_dict_1["ltc_param"], - [(list_modes[2]/list_modes[3]), (list_modes[2]/list_modes[3])]) - assert chromophore_dict_1["static_filter_list"] == ['LongEdge', [[True, False], kmax2]] - - # Case 2: list_lop, list_modes_by_bath, nmode_LTC, and static_filter_list - bath_dict_2 = {"list_lop": list_lop, "list_modes_by_bath": list_modes_by_bath, "nmodes_LTC": nmodes_LTC, "static_filter_list": static_filter_list_triangular} - chromophore_dict_2 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict_2) + assert np.allclose(chromophore_dict_1["ltc_param"], [0, 0]) + assert chromophore_dict_1["static_filter_list"] == [ + ['Markovian', [True, False, True, False]], + ['Triangular', [[True, False, True, False], kmax2]], + ['LongEdge', [[True, False, True, False], kmax2]]] + + # Case 2: list_lop, list_modes_by_bath, and static_filter_list + bath_dict_2 = {"list_lop": list_lop, "list_modes_by_bath": list_modes_by_bath, + "static_filter_list": static_filter_list_lmbb} + chromophore_dict_2 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + bath_dict_2) - # Note: when nmodes_LTC is set to 1, lop_list_hier=list_lop=lop_list_ltc for a dimer - for a, b in zip(chromophore_dict_2["lop_list_hier"], list_lop): + # Note: when nmodes_LTC is not set, lop_list_hier=lop_list_noise and lop_list_ltc=[] + for a, b in zip(chromophore_dict_2["lop_list_hier"], lop_list_noise): assert np.allclose(a.toarray(), b.toarray()) - for a, b in zip(chromophore_dict_2["lop_list_ltc"], list_lop): + for a, b in zip(chromophore_dict_2["lop_list_ltc"], []): assert np.allclose(a.toarray(), b.toarray()) - for a, b in zip(chromophore_dict_2["lop_list_noise"], lop_list_noise_nmodes_LTC): + for a, b in zip(chromophore_dict_2["lop_list_noise"], lop_list_noise): assert np.allclose(a.toarray(), b.toarray()) - assert np.allclose(chromophore_dict_2["gw_sysbath_hier"], [(list_modes_by_bath[0][0], - list_modes_by_bath[0][1]), - (list_modes_by_bath[1][0], - list_modes_by_bath[1][1])]) - assert np.allclose(chromophore_dict_2["gw_sysbath_noise"], [(list_modes_by_bath[0][0], - list_modes_by_bath[0][1]), - (list_modes_by_bath[0][2], - list_modes_by_bath[0][3]), - (list_modes_by_bath[1][0], - list_modes_by_bath[1][1]), - (list_modes_by_bath[1][2], - list_modes_by_bath[1][3])]) - assert np.allclose(chromophore_dict_2["ltc_param"], - [(list_modes_by_bath[0][2]/list_modes_by_bath[0][3]), - (list_modes_by_bath[1][2]/list_modes_by_bath[1][3])]) - assert chromophore_dict_2["static_filter_list"] == ['Triangular', [[True, False], kmax2]] - - # Case 3: list_modes only - bath_dict_3 = {"list_modes": list_modes} - chromophore_dict_3 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, bath_dict_3) - - for a, b in zip(chromophore_dict_3["lop_list_hier"], - [sparse.coo_matrix(([1], ([1], [1])), shape=(3, 3)), - sparse.coo_matrix(([1], ([1], [1])), shape=(3, 3)), - sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3)), - sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3))]): + assert np.allclose(chromophore_dict_2["gw_sysbath_hier"], + [(list_modes_by_bath[0][0], list_modes_by_bath[0][1]), + (list_modes_by_bath[0][2], list_modes_by_bath[0][3]), + (list_modes_by_bath[1][0], list_modes_by_bath[1][1]), + (list_modes_by_bath[1][2], list_modes_by_bath[1][3])]) + assert np.allclose(chromophore_dict_2["gw_sysbath_noise"], + [(list_modes_by_bath[0][0], list_modes_by_bath[0][1]), + (list_modes_by_bath[0][2], list_modes_by_bath[0][3]), + (list_modes_by_bath[1][0], list_modes_by_bath[1][1]), + (list_modes_by_bath[1][2], list_modes_by_bath[1][3])]) + assert np.allclose(chromophore_dict_2["ltc_param"],[0, 0]) + assert chromophore_dict_2["static_filter_list"] == [ + ['Markovian', [True, False, True, False]], + ['Triangular', [[True, False, True, False], kmax2]], + ['LongEdge', [[True, False, True, False], kmax2]]] + + # Case 3: list_modes + list_lop + nmodes_LTC + bath_dict_3 = {"list_modes": list_modes, "list_lop": list_lop, + "nmodes_LTC": nmodes_LTC} + chromophore_dict_3 = prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + bath_dict_3) + + for a, b in zip(chromophore_dict_3["lop_list_hier"], list_lop): assert np.allclose(a.toarray(), b.toarray()) - for a, b in zip(chromophore_dict_3["lop_list_noise"], - [sparse.coo_matrix(([1], ([1], [1])), shape=(3, 3)), - sparse.coo_matrix(([1], ([1], [1])), shape=(3, 3)), - sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3)), - sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3))]): + for a, b in zip(chromophore_dict_3["lop_list_noise"], lop_list_noise): assert np.allclose(a.toarray(), b.toarray()) - for a, b in zip(chromophore_dict_3["lop_list_ltc"], - [sparse.coo_matrix(([1], ([1], [1])), shape=(3, 3)), - sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3))]): + for a, b in zip(chromophore_dict_3["lop_list_ltc"], list_lop): assert np.allclose(a.toarray(), b.toarray()) assert np.allclose(chromophore_dict_3["gw_sysbath_hier"], [(list_modes[0], list_modes[1]), - (list_modes[2], - list_modes[3]), (list_modes[0], - list_modes[1]), - (list_modes[2], - list_modes[3])]) + list_modes[1])]) assert np.allclose(chromophore_dict_3["gw_sysbath_noise"], [(list_modes[0], list_modes[1]), (list_modes[2], @@ -742,174 +722,228 @@ def test_prepare_chromophore_input_dict(): list_modes[1]), (list_modes[2], list_modes[3])]) - assert np.allclose(chromophore_dict_3["ltc_param"], [0, 0]) + assert np.allclose(chromophore_dict_3["ltc_param"], [list_modes[2]/list_modes[3], + list_modes[2]/list_modes[3]]) assert chromophore_dict_3["static_filter_list"] == None - ## Testing M2_mu_ge input errors + # Testing M2_mu_ge input errors M2_mu_ge_wrongshape = np.array([np.array([0.5, 0.2]), np.array([0.5, 0.2])]) - try: - prepare_chromophore_input_dict(M2_mu_ge_wrongshape, H2_sys_hamiltonian, bath_dict_1) - except ValueError as excinfo: - if 'M2_mu_ge must be a numpy array with shape (n_chromophore, 3).' not in str(excinfo): - pytest.fail() + with pytest.raises(ValueError, match='M2_mu_ge must be a numpy array with shape ' + '\\(n_chromophore, 3\\).'): + prepare_chromophore_input_dict(M2_mu_ge_wrongshape, H2_sys_hamiltonian, + bath_dict_1) - ## Testing nmodes_LTC input errors + # Testing nmodes_LTC input errors nmodes_LTC_wrongtype = '1' nmodes_LTC_wrongvalue = -2 nmodes_LTC_large = 2 - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes,"nmodes_LTC": nmodes_LTC_wrongtype}) - except ValueError as excinfo: - if 'nmodes_LTC must be an integer or None.' not in str(excinfo): - pytest.fail() - - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "nmodes_LTC": nmodes_LTC_wrongvalue}) - except ValueError as excinfo: - if 'nmodes_LTC must be >= 0 or set to None.' not in str(excinfo): - pytest.fail() - - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes,"nmodes_LTC": nmodes_LTC_large}) - except ValueError as excinfo: - if 'nmodes_LTC must be less than the number of modes in each bath.' not in str(excinfo): - pytest.fail() - - ## Testing static_filter_list input errors + with pytest.raises(ValueError, match='nmodes_LTC must be an integer or None.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, { + "list_modes": list_modes,"nmodes_LTC": nmodes_LTC_wrongtype}) + + with pytest.raises(ValueError, match='nmodes_LTC must be >= 0 or set to None.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, { + "list_modes": list_modes, "nmodes_LTC": nmodes_LTC_wrongvalue}) + + with pytest.raises(ValueError, match='nmodes_LTC must be less than the number of ' + 'modes in each bath.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, { + "list_modes": list_modes,"nmodes_LTC": nmodes_LTC_large}) + + # Testing static_filter_list input errors # Case 1: static_filter_list not a list static_filter_list_wrongtype = 'Markovian' - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongtype}) - except ValueError as excinfo: - if 'static_filter_list must be a list.' not in str(excinfo): - pytest.fail() + with pytest.raises(ValueError, match='static_filter_list must be a list.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, { + "list_modes": list_modes, "static_filter_list": + static_filter_list_wrongtype}) # Case 2: static_filter_list not a list of length 2 [filter_name, filter_params] - static_filter_list_wronglength = ['Markovian'] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wronglength}) - except ValueError as excinfo: - if 'static_filter_list must be a 2-element list of the form: [filter_name, filter_params].' not in str(excinfo): - pytest.fail() - - # Case 3: static_filter_list[0] not an allowed option ["Markovian", "Triangular", "LongEdge"] - static_filter_list_wrongoption = ['Sharkovian', [True, False]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongoption}) - except ValueError as excinfo: - if "Filter name must be 'Markovian', 'Triangular', or 'LongEdge'." not in str(excinfo): - pytest.fail() - - # Case 4: filter_params not a list of length 2 [filter_bool, kmax2] for "Triangular" or "LongEdge" - static_filter_list_wrongparams = ['Triangular', [True, False, 4]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongparams}) - except ValueError as excinfo: - if 'The Triangular and LongEdge filter_params must be a list of booleans and an integer.' not in str(excinfo): - pytest.fail() + static_filter_list_wronglength = [['Markovian']] + with pytest.raises(ValueError, + match='static_filter_list must be a 2-element list of the form: ' + '\\[filter_name, filter_params\\].'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_wronglength}) + + # Case 3: static_filter_list[0] not an allowed option ["Markovian", + # "Triangular", "LongEdge"] + static_filter_list_wrongoption = [['Sharkovian', [True, False]]] + with pytest.raises(ValueError, + match="Error in filter 0: Filter names must be 'Markovian', " + "'Triangular', or 'LongEdge'."): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_wrongoption}) + + # Case 4: filter_params not a list of length 2 [filter_bool, kmax2] for + # "Triangular" or "LongEdge" + static_filter_list_wrongparams = [['Triangular', [True, False, 4]]] + with pytest.raises(ValueError, + match='Error in filter 0: Triangular/LongEdge filter_params ' + 'must be a list containing a list of booleans and an ' + 'integer.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_wrongparams}) # Case 5: Second entry of filter_params not an integer for "Triangular" or "LongEdge" - static_filter_list_wrongparams = ['Triangular', [[True, False], 'False']] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongparams}) - except ValueError as excinfo: - if 'The second entry in filter_params for the Triangular and LongEdge filters must be an integer.' not in str(excinfo): - pytest.fail() - - # Case 6: First entry of filter_params not a list of booleans for "Triangular" or "LongEdge" - static_filter_list_wrongparams = ['Triangular', [[True, 'False'], 2]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongparams}) - except ValueError as excinfo: - if 'The first entry in filter_params for the Triangular and LongEdge filters must be a list of booleans.' not in str(excinfo): - pytest.fail() - - # Case 7: kmax2 for "Triangular" or "LongEdge" not a positive integer - static_filter_list_wrongparams = ['Triangular', [[True, False], -2]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongparams}) - except ValueError as excinfo: - if 'Triangular and LongEdge filter_params must have a positive integer as the second element.' not in str(excinfo): - pytest.fail() + static_filter_list_wrongparams = [['Triangular', [[True, False], 'False']]] + with pytest.raises(ValueError, + match='Error in filter 0: The second entry in filter_params for ' + 'Triangular/LongEdge filters must be an integer.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": static_filter_list_wrongparams}) + + # Case 6: First entry of filter_params not a list of booleans for "Triangular" or + # "LongEdge" + static_filter_list_wrongparams = [['Triangular', [[True, 'False'], 2]]] + with pytest.raises(ValueError, + match='Error in filter 0: The first entry in filter_params for ' + 'Triangular/LongEdge filters must be a list of booleans.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_wrongparams}) + + # Case 7: Direct test for the nested conditional block - Triangular filter with + # negative kmax2 + static_filter_list_triangular_negative = [['Triangular', [[True, False], + kmax2_negative]]] + with pytest.raises(ValueError, + match='Error in filter 0: Triangular/LongEdge filter_params must ' + 'have a positive integer as the second element.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_triangular_negative}) # Case 8: static_filter_list not the right length for "Markovian" case - static_filter_list_wronglength = ['Markovian', [True]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wronglength}) - except ValueError as excinfo: - if 'The list of booleans in static_filter_list must have the same length as the number of modes.' not in str(excinfo): - pytest.fail() - - # Case 9: static_filter_list not the right length for "Triangular"/"LongEdge" cases - static_filter_list_wronglength = ['Triangular', [[True], 2]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wronglength}) - except ValueError as excinfo: - if 'The list of booleans in static_filter_list must have the same length as the number of modes.' not in str(excinfo): - pytest.fail() - - # Case 10: Filter_params not a list of booleans for "Markovian" - static_filter_list_wrongparams = ['Markovian', ['True', 'False']] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "static_filter_list": static_filter_list_wrongparams}) - except ValueError as excinfo: - if 'filter_params for Markovian filter must be a list of booleans.' not in str(excinfo): - pytest.fail() - - ## Testing that overdefining modes (list_modes_by_bath and list_modes) yields an error - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes, "list_modes_by_bath": list_modes_by_bath}) - except ValueError as excinfo: - if 'list_modes_by_bath and list_modes should not both be defined.' not in str(excinfo): - pytest.fail() - - ## Testing incompatible list_modes_by_bath and list_lop lengths yields an error + static_filter_list_wronglength = [['Markovian', [True]]] + with pytest.raises(ValueError, + match='Error in filter 0: The list of booleans in filter_params ' + 'must have the same length as the number of modes in each ' + 'bath.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_wronglength}) + + # Case 9: static_filter_list not the right length for "Markovian" case + # (list_modes_by_bath) + static_filter_list_wronglength = [['Markovian', [True, False]]] + with pytest.raises(ValueError, + match='Error in filter 0: The list of booleans in filter_params ' + 'must have the same length as the number of modes in all ' + 'baths combined.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes_by_bath": list_modes_by_bath, + "static_filter_list": + static_filter_list_wronglength}) + + # Case 10: static_filter_list not the right length for "Triangular"/"LongEdge" cases + static_filter_list_wronglength = [['Triangular', [[True], 2]]] + with pytest.raises(ValueError, + match='Error in filter 0: The list of booleans in filter_params ' + 'must have the same length as the number of modes in each ' + 'bath.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": + static_filter_list_wronglength}) + + # Case 11: static_filter_list not the right length for "Triangular"/"LongEdge" cases + # (list_modes_by_bath) + static_filter_list_wronglength = [['Triangular', [[True, False], 2]]] + with pytest.raises(ValueError, + match='Error in filter 0: The list of booleans in filter_params ' + 'must have the same length as the number of modes in all ' + 'baths combined.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes_by_bath": list_modes_by_bath, + "static_filter_list": + static_filter_list_wronglength}) + # Case 11: Filter_params not a list of booleans for "Markovian" + static_filter_list_wrongparams = [['Markovian', ['True', 'False']]] + with pytest.raises(ValueError, + match='Error in filter 0: filter_params for Markovian filters ' + 'must be a list of booleans.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "static_filter_list": static_filter_list_wrongparams}) + + # Testing that defining both nmodes_LTC and static_filter_list raises an error + with pytest.raises(ValueError, + match='The use of static hierarchy filters with low-temperature ' + 'correction is not currently supported.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "nmodes_LTC": nmodes_LTC, + "static_filter_list": static_filter_list_lm}) + + # Testing that overdefining modes (list_modes_by_bath and list_modes) yields error + with pytest.raises(ValueError, + match='list_modes_by_bath and list_modes should not both be ' + 'defined.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes": list_modes, + "list_modes_by_bath": list_modes_by_bath}) + + # Testing incompatible list_modes_by_bath and list_lop lengths yields an error list_modes_by_bath_wronglength = [ishizaki_decomposition_bcf_dl(40, 60, 290, 0)] list_lop_wronglength = [sparse.coo_matrix(([1], ([1], [2])), shape=(3, 3)), - sparse.coo_matrix(([1], ([2], [1])), shape=(3, 3)), sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3))] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_lop": list_lop_wronglength, "list_modes_by_bath": list_modes_by_bath_wronglength}) - except ValueError as excinfo: - if 'list_modes_by_bath and list_lop must have the same length.' not in str(excinfo): - pytest.fail() - - ## Testing list_modes_by_bath structure + sparse.coo_matrix(([1], ([2], [1])), shape=(3, 3)), + sparse.coo_matrix(([1], ([2], [2])), shape=(3, 3))] + with pytest.raises(ValueError, match='list_modes_by_bath and list_lop must have the' + ' same length.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_lop": list_lop_wronglength, + "list_modes_by_bath": + list_modes_by_bath_wronglength}) + + # Testing list_modes_by_bath structure # list_modes_by_bath not a list of lists list_modes_by_bath_wrongstructure = ['mode1', 'mode2'] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes_by_bath": list_modes_by_bath_wrongstructure}) - except ValueError as excinfo: - if 'list_modes_by_bath must be a list of lists.' not in str(excinfo): - pytest.fail() + with pytest.raises(ValueError, + match='list_modes_by_bath must be a list of lists.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes_by_bath": + list_modes_by_bath_wrongstructure}) # list_modes_by_bath doesn't contain paired Gs and Ws list_modes_by_bath_wrongpairing = [ishizaki_decomposition_bcf_dl(35, 50, 295, 0), [1, 2, 3]] - try: - prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes_by_bath": list_modes_by_bath_wrongpairing}) - except ValueError as excinfo: - if 'list_modes_by_bath should contain paired Gs and Ws, which guarantees an even number of elements in each sublist.' not in str(excinfo): - pytest.fail() - - ## Testing list_modes structure + with pytest.raises(ValueError, + match='list_modes_by_bath should contain paired Gs and Ws, which ' + 'guarantees an even number of elements in each sublist.'): + prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, + {"list_modes_by_bath": + list_modes_by_bath_wrongpairing}) + + # Testing list_modes structure # list_modes doesn't contain paired Gs and Ws list_modes_wrongpairing = [1, 2, 3] - try: + with pytest.raises(ValueError, match='list_modes should contain paired Gs and Ws, ' + 'which guarantees an even number of elements.'): prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {"list_modes": list_modes_wrongpairing}) - except ValueError as excinfo: - if 'list_modes should contain paired Gs and Ws, which guarantees an even number of elements.' not in str(excinfo): - pytest.fail() ## Testing that an error is raised if neither list_modes nor list_modes_by_bath are defined - try: + with pytest.raises(ValueError, match='Either list_modes_by_bath or list_modes must ' + 'be defined.'): prepare_chromophore_input_dict(M2_mu_ge, H2_sys_hamiltonian, {}) - except ValueError as excinfo: - if 'Either list_modes_by_bath or list_modes must be defined.' not in str(excinfo): - pytest.fail() def test_prepare_convergence_parameter_dict(): - ## Test proper output with all parameters defined + """ + Tests the prepare_convergence_parameter_dict helper function, ensuring that the + input dictionary is properly formatted with default and custom parameters. + """ + # Test proper output with all parameters defined t_step = 0.1 max_hier = 12 delta_a = 1e-10 @@ -927,11 +961,11 @@ def test_prepare_convergence_parameter_dict(): assert convergence_dict["set_update_step"] == set_update_step assert convergence_dict["set_f_discard"] == set_f_discard - ## Test proper output with defaults + # Test proper output with defaults convergence_dict = prepare_convergence_parameter_dict(t_step, max_hier) assert convergence_dict["t_step"] == t_step assert convergence_dict["max_hier"] == max_hier assert convergence_dict["delta_a"] == 0 assert convergence_dict["delta_s"] == 0 assert convergence_dict["set_update_step"] == 1 - assert convergence_dict["set_f_discard"] == 0 \ No newline at end of file + assert convergence_dict["set_f_discard"] == 0 diff --git a/tests/test_eom_hops_ksuper.py b/tests/test_eom_hops_ksuper.py index 1d42857..f77803b 100644 --- a/tests/test_eom_hops_ksuper.py +++ b/tests/test_eom_hops_ksuper.py @@ -11,7 +11,7 @@ ) from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp __title__ = "Test of eom_hops_ksuper" __author__ = "D. I. G. Bennett, B. Citty, J. K. Lynd" @@ -628,7 +628,7 @@ def test_add_crossterms_arbitrary_lop(): e_lambda = 65. gamma = 53.0 temp = 300.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) numsites2 = nsite*(nsite+1)/2 gw_sysbath = [] for i in range(nsite): @@ -914,7 +914,7 @@ def test_add_crossterms_stable_arbitrary_lop(): e_lambda = 65. gamma = 53.0 temp = 300.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) numsites2 = nsite*(nsite+1)/2 gw_sysbath = [] for i in range(nsite): diff --git a/tests/test_git_utils.py b/tests/test_git_utils.py new file mode 100644 index 0000000..e24a91e --- /dev/null +++ b/tests/test_git_utils.py @@ -0,0 +1,44 @@ +import pytest +from unittest.mock import patch, MagicMock +from mesohops.util.git_utils import get_git_commit_hash + + +def test_get_git_commit_hash_success(): + """Test that get_git_commit_hash returns a valid hash in a git repository.""" + # Since we're running in a git repository, we should get a valid hash + result = get_git_commit_hash() + + # Check that the result is a non-empty string (a valid git hash) + assert isinstance(result, str) + assert len(result) > 0 + + # Verify it's a valid git hash format (40 hex characters) + assert all(c in '0123456789abcdef' for c in result.lower()) + assert len(result) == 40 + + +@patch('subprocess.run') +def test_get_git_commit_hash_command_error(mock_run): + """Test that get_git_commit_hash returns error message when git command fails.""" + # Mock the subprocess.run to return a generic error + mock_process = MagicMock() + mock_process.returncode = 1 # Generic error code + mock_process.stderr = "A gigantic walrus ate the git repository" + mock_run.return_value = mock_process + + # Call the function and check the result + result = get_git_commit_hash() + assert result == "A gigantic walrus ate the git repository" + + +@patch('subprocess.run') +def test_get_git_commit_hash_file_not_found(mock_run): + """Test that get_git_commit_hash handles FileNotFoundError correctly.""" + # Mock subprocess.run to raise FileNotFoundError + mock_run.side_effect = FileNotFoundError + + # Call the function and check the result + result = get_git_commit_hash() + + # The function should return string when FileNotFoundError is caught + assert result == "Git is not installed or not found in PATH." diff --git a/tests/test_hierarchy_class.py b/tests/test_hierarchy_class.py index 4f6cae2..b36be59 100644 --- a/tests/test_hierarchy_class.py +++ b/tests/test_hierarchy_class.py @@ -1,12 +1,16 @@ +import pytest import numpy as np from mesohops.basis.hops_aux import AuxiliaryVector as AuxVec from mesohops.basis.hops_hierarchy import HopsHierarchy as HHier +from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS +from mesohops.trajectory.exp_noise import bcf_exp +from mesohops.util.exceptions import UnsupportedRequest __title__ = "test of hops hierarchy" __author__ = "D. I. G. Bennett, L. Varvelo, J. K. Lynd" -__version__ = "1.2" -__date__ = "Jan. 14, 2020" +__version__ = "1.6" +__date__ = "Aug. 13, 2025" def test_hierarchy_initialize_true(): @@ -324,6 +328,79 @@ def map_to_auxvec(list_aux, n_hmodes): ], 2 ) +def test_static_filter_init(): + """ + Tests that the static hierarchy filters throw an error if defined for the wrong + number of modes. + """ + # Define a HOPS trajectory in a non-special case + noise_param = {"SEED": None, "MODEL": "FFT_FILTER", "TLEN": 250.0, + # Units: fs + "TAU": 1.0, # Units: fs + } + nsite = 10 + (g_0, w_0) = [10,150] + loperator = np.zeros([10, 10, 10], dtype=np.float64) + gw_sysbath = [] + lop_list = [] + for i in range(nsite): + loperator[i, i, i] = 1.0 + gw_sysbath.append([g_0, w_0]) + lop_list.append(loperator[i]) + gw_sysbath.append([-1j * np.imag(g_0), 500.0]) + lop_list.append(loperator[i]) + hs = np.zeros([nsite, nsite]) + psi_0 = np.array([0.0] * nsite, dtype=np.complex128) + psi_0[5] = 1.0 + psi_0 = psi_0 / np.linalg.norm(psi_0) + + sys_param = {"HAMILTONIAN": np.array(hs, dtype=np.complex128), + "GW_SYSBATH": gw_sysbath, "L_HIER": lop_list, "L_NOISE1": lop_list, + "ALPHA_NOISE1": bcf_exp, "PARAM_NOISE1": gw_sysbath, } + + eom_param = {"EQUATION_OF_MOTION": "NORMALIZED NONLINEAR"} + + integrator_param = {"INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 5, + 'INCHWORM_CAP': 5, 'STATIC_BASIS': None} + + hierarchy_param_mark = {'MAXHIER': 6, + 'STATIC_FILTERS': [ + ['Markovian', [True] * (len(gw_sysbath) - 1)], + ]} + + hierarchy_param_tri = {'MAXHIER': 6, + 'STATIC_FILTERS': [ + ['Triangular', [[True] * (len(gw_sysbath) - 1), 2]], + ]} + + hierarchy_param_le = {'MAXHIER': 6, + 'STATIC_FILTERS': [ + ['LongEdge', [[True] * (len(gw_sysbath) - 1), 2]], + ]} + + for filter_dict in [hierarchy_param_mark, hierarchy_param_tri, hierarchy_param_le]: + with pytest.raises(UnsupportedRequest, match="The number of entries in the list " + "of static filter booleans does not"): + hops = HOPS(sys_param, + noise_param=noise_param, + hierarchy_param=filter_dict, + eom_param=eom_param, + integration_param=integrator_param, ) + hops.make_adaptive(0, 0) + hops.initialize(psi_0) + + with pytest.raises(UnsupportedRequest, match="The number of entries in the list " + "of static filter booleans does not"): + hops_ad = HOPS(sys_param, + noise_param=noise_param, + hierarchy_param=filter_dict, + eom_param=eom_param, + integration_param=integrator_param, ) + hops_ad.make_adaptive(1e-3, 1e-3) + hops_ad.initialize(psi_0) + hops_ad.propagate(10.0, 2.0) def test_define_triangular_hierarchy_4modes_4maxhier(): """ @@ -337,7 +414,6 @@ def test_define_triangular_hierarchy_4modes_4maxhier(): assert set(HH.define_triangular_hierarchy(4, 4)) == set(aux_list_4_4) assert set(HH.define_triangular_hierarchy(2, 4)) == set(aux_list_2_4) - def test_define_markovian_filtered_triangular_hierarchy(): """ Tests that the define_markovian_filtered_triangular_hierarchy function that is @@ -349,7 +425,7 @@ def test_define_markovian_filtered_triangular_hierarchy(): hierarchy_param = {"MAXHIER": 10, "STATIC_FILTERS": [['Markovian', [False, True, True, True]]]} - system_param = {"N_HMODES": 10} + system_param = {"N_HMODES": 4} HH = HHier(hierarchy_param, system_param) mark_list_aux = HH.define_markovian_filtered_triangular_hierarchy(4, 4, [False, True, True, True]) diff --git a/tests/test_hops_aux.py b/tests/test_hops_aux.py index 3c1ef70..15f7350 100644 --- a/tests/test_hops_aux.py +++ b/tests/test_hops_aux.py @@ -1,4 +1,6 @@ import numpy as np +import pytest + from mesohops.basis.hops_aux import AuxiliaryVector from mesohops.util.exceptions import AuxError @@ -10,10 +12,8 @@ def test_auxvec_ordering(): """ aux_1010 = AuxiliaryVector([(0, 1), (2, 1)], 4) assert type(aux_1010) == AuxiliaryVector - try: + with pytest.raises(AuxError, match='array_aux_vec not properly ordered'): aux_1010 = AuxiliaryVector([(2, 1), (0, 1)], 4) - except AuxError as excinfo: - assert 'array_aux_vec not properly ordered' in str(excinfo) def test_keys(): @@ -221,11 +221,9 @@ def test_add_aux_connect(): assert aux_1001.dict_aux_m1[3] == aux_1000 # # Test when type != +/- 1 - try: + with pytest.raises(AuxError, match='There is a problem in the hierarchy: ' + 'add_aux_connect does not support type=2'): aux_1002.add_aux_connect(3, aux_1000, 2) - except AuxError as excinfo: - assert 'There is a problem in the hierarchy: add_aux_connect does not support ' \ - 'type=2' in str(excinfo) def test_remove_aux_connect(): @@ -250,12 +248,9 @@ def test_remove_aux_connect(): assert aux_1001.dict_aux_m1 == {} # Test when type != +/- 1 - try: + with pytest.raises(AuxError, match='There is a problem in the hierarchy: ' + 'remove_aux_connect does not support type=2'): aux_1001.remove_aux_connect(3, 2) - except AuxError as excinfo: - assert 'There is a problem in the hierarchy: remove_aux_connect does not ' \ - 'support ' \ - 'type=2' in str(excinfo) def test_remove_pointers(): """ diff --git a/tests/test_hops_basis.py b/tests/test_hops_basis.py index 570774a..8142b07 100644 --- a/tests/test_hops_basis.py +++ b/tests/test_hops_basis.py @@ -1,12 +1,14 @@ -import pytest import os + import numpy as np +import pytest import scipy as sp + from mesohops.basis.hops_aux import AuxiliaryVector as AuxiliaryVector from mesohops.basis.hops_hierarchy import HopsHierarchy as HHier from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.util.exceptions import UnsupportedRequest from mesohops.util.physical_constants import hbar @@ -57,7 +59,7 @@ def test_initialize(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -172,7 +174,7 @@ def test_define_basis_state(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -268,7 +270,7 @@ def test_define_basis_hier(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([2, 2, 2], dtype=np.float64) gw_sysbath = [] @@ -362,7 +364,7 @@ def test_update_basis(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -490,7 +492,7 @@ def test_define_state_basis(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -576,11 +578,10 @@ def test_define_state_basis(): known_boundary = [3, 7] assert np.array_equal(list_state_bound, known_boundary) - try: - list_stable_state, list_state_bound = hops_ad.basis._define_state_basis( - phi_new/5, 2.0, z_step, list_index_aux_stable, []) - except UnsupportedRequest as excinfo: - assert "does not support non-normalized Φ in the _define_state" in str(excinfo) + with pytest.raises(UnsupportedRequest, match="does not support non-normalized Φ in " + "the _define_state"): + hops_ad.basis._define_state_basis(phi_new/5, 2.0, z_step, + list_index_aux_stable, []) def test_define_hierarchy_basis(): """ @@ -598,7 +599,7 @@ def test_define_hierarchy_basis(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([2, 2, 2], dtype=np.float64) gw_sysbath = [] @@ -679,13 +680,11 @@ def test_define_hierarchy_basis(): assert set(list_aux_boundary) == set(known_boundary) - try: + with pytest.raises(UnsupportedRequest, match="does not support non-normalized Φ in " + "the _define_hierarchy"): list_aux_stable, list_aux_boundary = hops_ad.basis._define_hierarchy_basis( phi_test/5, 2.0, z_step ) - except UnsupportedRequest as excinfo: - assert "does not support non-normalized Φ in the _define_hierarchy" in str( - excinfo) def test_determine_boundary_hier(): @@ -703,7 +702,7 @@ def test_determine_boundary_hier(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([2, 2, 2], dtype=np.float64) gw_sysbath = [] @@ -995,7 +994,7 @@ def test_fraction_discard(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([2, 2, 2], dtype=np.float64) gw_sysbath = [] @@ -1041,16 +1040,13 @@ def test_fraction_discard(): eom_param=eom_param, integration_param=integrator_param, ) - try: + with pytest.raises(UnsupportedRequest) as excinfo: hops_ad.make_adaptive(1e-3, 1e-3, f_discard=-1) - except UnsupportedRequest as excinfo: - if "acceptable range of [0,1]" not in str(excinfo): - pytest.fail() - try: + assert "acceptable range of [0,1]" in str(excinfo.value) + + with pytest.raises(UnsupportedRequest) as excinfo: hops_ad.make_adaptive(1e-3, 1e-3, f_discard=2) - except UnsupportedRequest as excinfo: - if "acceptable range of [0,1]" not in str(excinfo): - pytest.fail() + assert "acceptable range of [0,1]" in str(excinfo.value) # test that f_discard is properly assigned hops_ad = HOPS( @@ -1138,7 +1134,7 @@ def test_determine_basis_from_list(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -1242,7 +1238,7 @@ def test_state_stable_error(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -1370,7 +1366,7 @@ def test_list_M2_by_dest(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) gw_sysbath = [] def get_holstein(n): @@ -1661,7 +1657,7 @@ def test_get_Z2_noise_sparse(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) gw_sysbath = [] @@ -1771,7 +1767,7 @@ def get_peierls(n): Z2_noise_sparse_known = np.sum((np.array([noise_mem[m]*lop_list[m] for m in list_mode_off_diag])), axis=0) + np.sum(np.array( - [(np.conj(noise_1)-2j*np.real(noise_2))[m] * lop_list_base[m] for m in + [(np.conj(noise_1)-1j*noise_2)[m] * lop_list_base[m] for m in list_lop_in_basis_off_diag]), axis=0) Z2_noise_sparse = hops_ad.basis.get_Z2_noise_sparse(z_step) @@ -1813,7 +1809,7 @@ def test_get_T2_ltc(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) gw_sysbath = [] @@ -2000,7 +1996,7 @@ def test_hier_max_adaptive(): gamma = 60 temp = 300 hs = np.zeros([1, 1], np.complex128) - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) gw_sysbath = [[g_0, w_0]] lop_list = np.ones([1, 1, 1]) diff --git a/tests/test_hops_dyadic.py b/tests/test_hops_dyadic.py index bb00fce..6983d84 100644 --- a/tests/test_hops_dyadic.py +++ b/tests/test_hops_dyadic.py @@ -1,10 +1,12 @@ -import pytest -import numpy as np import time as timer + +import numpy as np +import pytest +from scipy import sparse + from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_dyadic import DyadicTrajectory as DHOPS from mesohops.util.exceptions import UnsupportedRequest -from scipy import sparse noise_param = { "SEED": 10, @@ -201,6 +203,20 @@ def test_M2_dyad_conversion(): np.testing.assert_allclose(sys_param_dense['L_NOISE1'], list(list_lop_dense_ref)*2) np.testing.assert_allclose(sys_param_dense['L_LT_CORR'], list_lop_dense_ref) + + M2_complex = np.array([[1 + 2j, 2 - 1j], + [3 + 0j, 4 - 0.5j]], dtype=np.complex128) + M2_dyad_ref = np.array([[1 + 2j, 2 - 1j, 0, 0], + [3 + 0j, 4 - 0.5j, 0, 0],[0, 0, 1 + 2j, 2 - 1j], + [0, 0, 3 + 0j, 4 - 0.5j]], dtype=np.complex128) + + M2_dyad_test = dhops_dense._M2_dyad_conversion(M2_complex) # call method manually + + np.testing.assert_allclose(M2_dyad_ref, M2_dyad_test) + # Check dtype + assert sys_param_dense['HAMILTONIAN'].dtype == np.complex128 + assert M2_dyad_test.dtype == np.complex128 + # Sparse Construct # ---------------- @@ -239,11 +255,9 @@ def test_dyad_operator(): Op_bra_dense @ psi_b))) ** 2)) np.testing.assert_allclose(dhops_dense.psi, psi_2_ref) - try: + with pytest.raises(UnsupportedRequest) as excinfo: dhops_dense._dyad_operator(Op_bra_dense, 'braket') - except UnsupportedRequest as excinfo: - if 'sides other than "ket" or "bra"' not in str(excinfo): - pytest.fail() + assert 'sides other than "ket" or "bra"' in str(excinfo.value) # Sparse Construct # ---------------- @@ -262,11 +276,9 @@ def test_dyad_operator(): Op_bra_sparse@psi_b))) ** 2)) np.testing.assert_allclose(dhops_sparse.psi, psi_2_ref) - try: + with pytest.raises(UnsupportedRequest) as excinfo: dhops_sparse._dyad_operator(Op_bra_sparse, 'braket') - except UnsupportedRequest as excinfo: - if 'sides other than "ket" or "bra"' not in str(excinfo): - pytest.fail() + assert 'sides other than "ket" or "bra"' in str(excinfo.value) @pytest.mark.order(4) def test_norm_comp_list(): diff --git a/tests/test_hops_eom.py b/tests/test_hops_eom.py index 01090eb..ba5392a 100644 --- a/tests/test_hops_eom.py +++ b/tests/test_hops_eom.py @@ -352,7 +352,7 @@ def test_linear_eom(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] list_noise_memory = hops.z_mem[hops.basis.mode.list_absindex_mode] @@ -379,7 +379,7 @@ def test_linear_eom(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] @@ -422,7 +422,7 @@ def test_normalized_nonlinear_nonadaptive_eom(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] list_noise_memory = hops.z_mem[hops.basis.mode.list_absindex_mode] @@ -449,7 +449,7 @@ def test_normalized_nonlinear_nonadaptive_eom(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] @@ -491,7 +491,7 @@ def test_nonlinear_eom(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] list_noise_memory = hops.z_mem[hops.basis.mode.list_absindex_mode] @@ -517,7 +517,7 @@ def test_nonlinear_eom(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] @@ -617,7 +617,7 @@ def test_eom_adaptive(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] @@ -648,7 +648,7 @@ def test_eom_adaptive(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [hops.basis.mode.list_g, hops.basis.mode.list_w, hops.basis.mode.list_absindex_mode] list_l_op = [sys_param["L_HIER"][m] for m in hops.basis.mode.list_absindex_mode] @@ -690,7 +690,7 @@ def test_eom_adaptive(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [linear_chain_hops.basis.mode.list_g, linear_chain_hops.basis.mode.list_w, linear_chain_hops.basis.mode.list_absindex_mode] list_l_op = [linear_chain_sys_param["L_HIER"][m] for m in @@ -726,7 +726,7 @@ def test_eom_adaptive(): noise2_t = np.arange(len(noise_t)) # Note: the noise we construct here is pre-conjugation, so noise 2 needs a # sign of +2j, rather than -2j. - noise_t_combined = noise_t + 2.0j * np.real(noise2_t) + noise_t_combined = noise_t + 1.0j * noise2_t list_modes = [linear_chain_hops.basis.mode.list_g, linear_chain_hops.basis.mode.list_w, linear_chain_hops.basis.mode.list_absindex_mode] list_l_op = [linear_chain_sys_param["L_HIER"][m] for m in diff --git a/tests/test_hops_fluxfilters.py b/tests/test_hops_fluxfilters.py index 4e2b2b1..79a4be7 100644 --- a/tests/test_hops_fluxfilters.py +++ b/tests/test_hops_fluxfilters.py @@ -1,7 +1,7 @@ import numpy as np import scipy as sp from mesohops.basis.hops_aux import AuxiliaryVector as AuxiliaryVector -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS @@ -15,7 +15,7 @@ def test_filter_hierarchy_stable_up(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -81,7 +81,7 @@ def test_filter_hierarchy_stable_down(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -147,7 +147,7 @@ def test_filter_hierarchy_boundary_up(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -211,7 +211,7 @@ def test_filter_hierarchy_boundary_down(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -274,7 +274,7 @@ def test_filter_state_stable_up(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -364,7 +364,7 @@ def test_filter_state_stable_down(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -442,7 +442,7 @@ def test_filter_markovian_up(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -523,6 +523,26 @@ def test_filter_markovian_up(): assert np.all(filter_markovian == known_filter_markovian) + # Test empty mode list + sys_param = {"HAMILTONIAN": np.array(hs, dtype=np.complex128), + "GW_SYSBATH": gw_sysbath[14:], + "L_HIER": lop_list[14:], + "L_NOISE1": lop_list[14:], + "ALPHA_NOISE1": bcf_exp, + "PARAM_NOISE1": gw_sysbath[14:], } + # Adaptive Hops + hops_ad = HOPS(sys_param, + noise_param=noise_param, + hierarchy_param={'MAXHIER': 6, + 'STATIC_FILTERS': [['Markovian', + list_mark_filter[14:]], + ['Markovian', list_mark_filter_2[14:]]], }, + eom_param=eom_param, + integration_param=integrator_param, ) + hops_ad.make_adaptive(1e-3, 1e-3) + hops_ad.initialize(psi_0) + assert hops_ad.basis.flux_filters.construct_filter_markov_up() == True + def test_filter_triangular_up(): """ Tests that the triangular static filter that filters subset of modes M with @@ -541,7 +561,7 @@ def test_filter_triangular_up(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -625,6 +645,26 @@ def test_filter_triangular_up(): assert np.all(filter_triangular == known_filter_triangular) + # Test empty mode list + sys_param = {"HAMILTONIAN": np.array(hs, dtype=np.complex128), + "GW_SYSBATH": gw_sysbath[14:], + "L_HIER": lop_list[14:], + "L_NOISE1": lop_list[14:], + "ALPHA_NOISE1": bcf_exp, + "PARAM_NOISE1": gw_sysbath[14:], } + # Adaptive Hops + hops_ad = HOPS(sys_param, + noise_param=noise_param, + hierarchy_param={'MAXHIER': 6, + 'STATIC_FILTERS': [ + ['Triangular', [list_tri_filter[14:], 2]], + ['Triangular', [list_tri_filter[14:], 3]]], }, + eom_param=eom_param, + integration_param=integrator_param, ) + hops_ad.make_adaptive(1e-3, 1e-3) + hops_ad.initialize(psi_0) + assert hops_ad.basis.flux_filters.construct_filter_triangular_up() == True + def test_filter_longedge_up(): """ Tests that the longedge static filter that filters subset of modes M with @@ -646,7 +686,7 @@ def test_filter_longedge_up(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -734,6 +774,26 @@ def test_filter_longedge_up(): assert np.all(filter_longedge == known_filter_longedge) + # Test empty mode list + sys_param = {"HAMILTONIAN": np.array(hs, dtype=np.complex128), + "GW_SYSBATH": gw_sysbath[14:], + "L_HIER": lop_list[14:], + "L_NOISE1": lop_list[14:], + "ALPHA_NOISE1": bcf_exp, + "PARAM_NOISE1": gw_sysbath[14:], } + # Adaptive Hops + hops_ad = HOPS(sys_param, + noise_param=noise_param, + hierarchy_param={'MAXHIER': 6, + 'STATIC_FILTERS': [ + ['LongEdge', [list_le_filter[14:], 2]], + ['LongEdge', [list_le_filter_2[14:], 3]]], }, + eom_param=eom_param, + integration_param=integrator_param, ) + hops_ad.make_adaptive(1e-3, 1e-3) + hops_ad.initialize(psi_0) + assert hops_ad.basis.flux_filters.construct_filter_longedge_up() == True + def test_mode_setter(): """ Tests that the setter for HopsModes.list_absindex_modes updates the relevant @@ -747,7 +807,7 @@ def test_mode_setter(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] diff --git a/tests/test_hops_mode.py b/tests/test_hops_mode.py index 1087d52..5cdd34d 100644 --- a/tests/test_hops_mode.py +++ b/tests/test_hops_mode.py @@ -3,7 +3,7 @@ from mesohops.basis.hops_aux import AuxiliaryVector as AuxiliaryVector from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp def test_mode_setter(): @@ -15,7 +15,7 @@ def test_mode_setter(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -112,7 +112,7 @@ def test_empty_modelist(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -169,7 +169,7 @@ def test_list_off_diag_active_mask(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) gw_sysbath = [] def get_holstein(n): diff --git a/tests/test_hops_noise.py b/tests/test_hops_noise.py index f5e677e..4f7901f 100644 --- a/tests/test_hops_noise.py +++ b/tests/test_hops_noise.py @@ -1,7 +1,8 @@ -import pytest import numpy as np -from mesohops.trajectory.exp_noise import bcf_exp +import pytest + from mesohops.noise.hops_noise import HopsNoise +from mesohops.trajectory.exp_noise import bcf_exp from mesohops.util.exceptions import UnsupportedRequest __title__ = "Test of hops_noise" @@ -94,11 +95,13 @@ def test_initialize(): assert noise_corr_working['N_L2'] == test_noise.param['N_L2'] assert noise_corr_working['LIND_BY_NMODE'] == test_noise.param['LIND_BY_NMODE'] assert noise_corr_working['CORR_PARAM'] == test_noise.param['CORR_PARAM'] + # Test that FLAG_REAL defaults to False + assert test_noise.param["FLAG_REAL"] == False # Test that the keys overlap excepting T_AXIS (added by HopsNoise) and # STORE_RAW_NOISE (added by FFTFilterNoise) assert set(list(noise_param.keys()) + list(noise_corr_working.keys()) + [ - 'T_AXIS', 'RAND_MODEL', 'STORE_RAW_NOISE', 'NOISE_WINDOW', 'ADAPTIVE']) == set( - test_noise.param.keys()) + 'T_AXIS', 'RAND_MODEL', 'STORE_RAW_NOISE', 'NOISE_WINDOW', 'ADAPTIVE', + 'FLAG_REAL' ]) == set(test_noise.param.keys()) def test_get_noise(capsys): @@ -162,12 +165,10 @@ def test_get_noise(capsys): # Tests that get_noise raises an UnsupportedRequest if using a nonexistent # NOISE_MODEL - try: + with pytest.raises(UnsupportedRequest) as excinfo: HopsNoise(noise_param_broken, noise_corr_working).get_noise(t_axis[:2]) - except UnsupportedRequest as excinfo: - if 'does not support Noise.param[MODEL] NONEXISTENT_NOISE in the ' not in str( - excinfo): - pytest.fail() + assert ('does not support Noise.param[MODEL] NONEXISTENT_NOISE in the ' in + str(excinfo.value)) # Windowed noise test test_noise_windowed = HopsNoise(noise_param_windowed, noise_corr_working) @@ -240,6 +241,58 @@ def test_get_noise(capsys): assert ("Warning: noise windowing is not supported while using interpolated " "noise") in out + # Tests of FLAG_REAL + noise_param_real = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 1000.0, # Units: fs + "TAU": 1.0, # Units: fs, + "INTERPOLATE": False, + "FLAG_REAL": True, + } + test_noise_real = HopsNoise(noise_param_real, noise_corr_working) + # if FLAG_REAL, noise should be purely real. + assert np.allclose(test_noise_real.get_noise(t_axis[:2])[0, :], + np.real(test_noise_real._noise[0, :2]), atol=1e-8) + + noise_param_complex = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 1000.0, # Units: fs + "TAU": 1.0, # Units: fs, + "INTERPOLATE": False, + "FLAG_REAL": False, + } + test_noise_complex = HopsNoise(noise_param_complex, noise_corr_working) + # If not FLAG_REAL, noise is not set to real. + assert not np.allclose(test_noise_complex.get_noise(t_axis[:2])[0, :], + np.real(test_noise_complex._noise[0, :2]), atol=1e-8) + + # Same test for interpolated noise + noise_param_interp_real = { + "SEED": 1j*noise, + "MODEL": "PRE_CALCULATED", + "TLEN": 1000.0, # Units: fs + "TAU": 1.0, # Units: fs, + "INTERPOLATE": True, + "FLAG_REAL": True, + } + test_noise_interp_real = HopsNoise(noise_param_interp_real, noise_corr_working) + assert np.allclose(test_noise_interp_real.get_noise([0, 0.25, 0.5, 0.75, 1]), + np.zeros([2,5]), atol=1e-8) + + noise_param_interp_complex = { + "SEED": 1j * noise, + "MODEL": "PRE_CALCULATED", + "TLEN": 1000.0, # Units: fs + "TAU": 1.0, # Units: fs, + "INTERPOLATE": True, + "FLAG_REAL": False, + } + test_noise_interp_complex = HopsNoise(noise_param_interp_complex, noise_corr_working) + assert not np.allclose(test_noise_interp_complex.get_noise([0, 0.25, 0.5, 0.75, 1]), + np.zeros([2, 5]), atol=1e-8) + def test_noise_adaptivity(): """ diff --git a/tests/test_hops_storage.py b/tests/test_hops_storage.py index 765c29b..6f0b721 100644 --- a/tests/test_hops_storage.py +++ b/tests/test_hops_storage.py @@ -1,10 +1,13 @@ import numpy as np +import pytest + import mesohops.storage.storage_functions as sf from mesohops.basis.hops_aux import AuxiliaryVector -from mesohops.trajectory.exp_noise import bcf_exp from mesohops.storage.hops_storage import HopsStorage +from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory from mesohops.util.exceptions import UnsupportedRequest +from unittest.mock import patch # New Hops_Storage tests @@ -35,10 +38,9 @@ def fake_save_func(): assert broken_HS.dic_save["t_axis"] == fake_save_func assert not "t_axis" in false_AHS.dic_save - try: - empty_AHS = HopsStorage(True, {"t_axis": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"t_axis": None}) def test_psi_traj_save(): @@ -61,10 +63,9 @@ def fake_save_func(): assert broken_HS.dic_save["psi_traj"] == fake_save_func assert not "psi_traj" in false_AHS.dic_save - try: - empty_AHS = HopsStorage(True, {"psi_traj": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"psi_traj": None}) def test_phi_traj_save(): @@ -85,10 +86,9 @@ def fake_save_func(): assert broken_HS.dic_save["phi_traj"] == fake_save_func assert not "phi_traj" in false_AHS.dic_save - try: - empty_AHS = HopsStorage(True, {"phi_traj": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"phi_traj": None}) def test_aux_new_save(): @@ -110,10 +110,9 @@ def fake_save_func(): assert broken_HS.dic_save["aux_list"] == fake_save_func assert not "aux_list" in false_AHS.dic_save - try: - empty_AHS = HopsStorage(True, {"aux_list": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"aux_list": None}) def test_state_list_save(): @@ -134,10 +133,9 @@ def fake_save_func(): assert broken_HS.dic_save["state_list"] == fake_save_func assert not "state_list" in false_AHS.dic_save.keys() - try: - empty_AHS = HopsStorage(True, {"state_list": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"state_list": None}) def test_list_nhier_save(): @@ -158,10 +156,9 @@ def fake_save_func(): assert broken_HS.dic_save["list_nhier"] == fake_save_func assert not "list_nhier" in false_AHS.dic_save.keys() - try: - empty_AHS = HopsStorage(True, {"list_nhier": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"list_nhier": None}) def test_list_nstate_save(): @@ -183,10 +180,9 @@ def fake_save_func(): assert broken_HS.dic_save["list_nstate"] == fake_save_func assert not "list_nstate" in false_AHS.dic_save.keys() - try: - empty_AHS = HopsStorage(True, {"list_nstate": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"list_nstate": None}) def test_list_aux_norm_save(): """ @@ -207,10 +203,9 @@ def fake_save_func(): assert broken_HS.dic_save["list_aux_norm"] == fake_save_func assert not "list_aux_norm" in false_AHS.dic_save.keys() - try: - empty_AHS = HopsStorage(True, {"list_aux_norm": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"list_aux_norm": None}) test_phi = [1,2,3,4,5,6,7,8] fake_aux_list = ['k0', 'k1', 'k2', 'k3'] @@ -237,10 +232,9 @@ def fake_save_func(): assert broken_HS.dic_save["z_mem"] == fake_save_func assert not "z_mem" in false_AHS.dic_save.keys() - try: - empty_AHS = HopsStorage(True, {"z_mem": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"z_mem": None}) fake_zmem_list = np.array([15+1j, 1234j-2, 234j-100, 137.23-0.2j]) test_phi = [1, 2, 3, 4, 5, 6, 7, 8] @@ -263,10 +257,9 @@ def dummy_saving_function(): assert HS.dic_save["property"] == dummy_saving_function assert AHS.dic_save["property"] == dummy_saving_function assert not "property" in false_AHS.dic_save - try: - empty_AHS = HopsStorage(True, {"property": None}) - except UnsupportedRequest as excinfo: - assert 'The current code does not support this value' in str(excinfo) + with pytest.raises(UnsupportedRequest, match='The current code does not support ' + 'this value'): + HopsStorage(True, {"property": None}) def test_store_step(): @@ -395,3 +388,137 @@ def test_hierarchy_storage_functions(): assert storage["list_nhier"][-1] == 4 for i in range(4): assert np.allclose(list_aux[i].array_aux_vec, storage["aux_list"][-1][i]) + +def test_check_storage_time(): + """ + Tests that the helper function that checks whether the data associated with a + given time point ought to be saved is working correctly. + """ + # Check that every time point is stored if STORAGE_TIME key not included. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True} + ) + assert hops_traj.storage.check_storage_time(-111) + + # Check that every time point is stored if STORAGE_TIME is True. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True, "STORAGE_TIME": True} + ) + assert hops_traj.storage.check_storage_time(-111) + + # Check that every time point is dropped if STORAGE_TIME is False. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True, "STORAGE_TIME": False} + ) + assert not hops_traj.storage.check_storage_time(-111) + + # Check that every time point evenly divisible by a integer STORAGE_TIME is saved. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True, "STORAGE_TIME": 10} + ) + assert hops_traj.storage.check_storage_time(10) + assert hops_traj.storage.check_storage_time(20) + assert not hops_traj.storage.check_storage_time(5) + + # Check that every time point evenly divisible by a float STORAGE_TIME is saved. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True, "STORAGE_TIME": 10.0} + ) + assert hops_traj.storage.check_storage_time(10) + assert hops_traj.storage.check_storage_time(20) + assert not hops_traj.storage.check_storage_time(5) + + # Check that every time point in an iterable STORAGE_TIME is saved. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True, "STORAGE_TIME": [5,10,25]} + ) + assert hops_traj.storage.check_storage_time(10) + assert not hops_traj.storage.check_storage_time(20) + assert hops_traj.storage.check_storage_time(5) + + # Integrated test that tests propagation has no indexing mistakes and always + # saves the 0th time point. Here, we should only save the 10 fs point, because 5 + # and 25 fs will never occur with a dt of 2.0 fs. + hops_traj.initialize(psi_0) + hops_traj.propagate(30.0, 2.0) + assert len(hops_traj.storage['list_nhier']) == 2 + + # Check that an error is raised when STORAGE_TIME is an unacceptable type. + hops_traj = HopsTrajectory( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param, + storage_param={"list_nhier": True, "aux_list": True, + "STORAGE_TIME": 'a'} + ) + with pytest.raises(UnsupportedRequest, match='is not allowed in the STORAGE_TIME ' + 'key of storage'): + hops_traj.storage.check_storage_time(10) + + +@patch('mesohops.storage.hops_storage.get_git_commit_hash') +def test_git_commit_hash_in_metadata(mock_get_hash): + """ + Test that the git commit hash is correctly stored in the metadata when a + HopsStorage object is created. + """ + # Mock the get_git_commit_hash function to return a test value + mock_get_hash.return_value = "No walruses were harmed in the making of this test" + + # Create a HopsStorage object + storage = HopsStorage(adaptive=False, storage_dic={}) + + # Check that the git commit hash is in the metadata + assert "GIT_COMMIT_HASH" in storage.metadata + assert storage.metadata["GIT_COMMIT_HASH"] == ("No walruses were harmed in the " + "making of this test") + + +@patch('mesohops.storage.hops_storage.get_git_commit_hash') +def test_git_commit_hash_error_in_metadata(mock_get_hash): + """ + Test that git error message is correctly stored in the metadata when git commit + hash cannot be retrieved. + """ + # Mock the get_git_commit_hash function to return an error message + mock_get_hash.return_value = "A gigantic walrus ate the git repository here too" + + # Create a HopsStorage object + storage = HopsStorage(adaptive=False, storage_dic={}) + + # Check that the git commit hash is in the metadata and is correct + assert "GIT_COMMIT_HASH" in storage.metadata + assert storage.metadata["GIT_COMMIT_HASH"] == ("A gigantic walrus ate the git " + "repository here too") diff --git a/tests/test_hops_system.py b/tests/test_hops_system.py index e13d57c..8ca81ba 100644 --- a/tests/test_hops_system.py +++ b/tests/test_hops_system.py @@ -1,9 +1,12 @@ import os +import pytest import numpy as np import scipy as sp from mesohops.basis.hops_system import HopsSystem as HSystem +from mesohops.basis.system_functions import initialize_system_dict from mesohops.trajectory.exp_noise import bcf_exp -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp +from .utils import compare_dictionaries __title__ = "test for System Class" __author__ = "D. I. G. Bennett, L. Varvelo" @@ -22,7 +25,7 @@ e_lambda = 20.0 gamma = 50.0 temp = 140.0 -(g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) +(g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([4, 4, 4], dtype=np.float64) gw_sysbath = [] @@ -139,28 +142,6 @@ def test_initialize_false(): assert np.array_equal(state, known_state) -def test_get_states_from_L2_try(): - """ - This function test to make sure _get_states_from_L2 is properly listing the states - that the L operators interacts with - """ - lop = lop_list[2] - lop = HS._get_states_from_L2(lop) - known_state_tuple = tuple([1]) - assert lop == known_state_tuple - - -def test_get_states_from_L2_except(): - """ - This function test to make sure _get_states_from_L2 is properly listing the states - that the L operators interacts with - """ - lop = sp.sparse.coo_matrix(lop_list[2]) - lop = HS._get_states_from_L2(lop) - known_state_tuple = tuple([1]) - assert lop == known_state_tuple - - def test_state_list_setter(): """ Tests that the state list setter correctly manages helper objects for indexing @@ -189,7 +170,7 @@ def test_state_list_setter(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([nsite, nsite, nsite], dtype=np.float64) gw_sysbath = [] lop_list = [] @@ -258,7 +239,7 @@ def test_state_list_setter(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) gw_sysbath = [] lop_list = [] for loperatori in list_loperator: @@ -394,3 +375,43 @@ def test_dict_relative_index_by_state(): assert HS.dict_relative_index_by_state[2] == 1 assert HS.dict_relative_index_by_state[3] == 2 + +def test_from_file_missing(): + """Ensures that constructing with a missing file raises a FileNotFoundError.""" + with pytest.raises(FileNotFoundError): + hs_loaded = HSystem("missing.pkl") + + +def test_save_and_load_full(tmp_path): + """Comprehensively tests saving and loading of system parameters.""" + # Small two state system used to keep the test light weight + ham = np.zeros((2, 2)) + l1 = np.array([[1.0, 0.0], [0.0, 0.0]]) + l2 = np.array([[0.0, 0.0], [0.0, 1.0]]) + l3 = np.array([[0.0, 1.0], [1.0, 0.0]]) + + param = { + "HAMILTONIAN": ham, + "GW_SYSBATH": [(0.1, 1.0), (0.2, 1.1), (0.3, 1.2)], + "L_HIER": [l1, l2, l3], + "L_NOISE1": [l1, l2, l3], + "PARAM_NOISE1": [(0.1, 1.0), (0.2, 1.1), (0.3, 1.2)], + "ALPHA_NOISE1": bcf_exp, + } + + hs = HSystem(param) + + fname = tmp_path / "sys.pkl" + hs.save_dict_param(fname) + assert fname.exists() + + hs_loaded = HSystem(fname) + + # Ensure all keys are present and values are equal + compare_dictionaries(hs.param, hs_loaded.param) + + +def test_load_invalid_type(): + """Ensures that constructing with an invalid type raises a TypeError.""" + with pytest.raises(TypeError): + HSystem(123) diff --git a/tests/test_hops_trajectory.py b/tests/test_hops_trajectory.py index cc2332b..3290e9e 100644 --- a/tests/test_hops_trajectory.py +++ b/tests/test_hops_trajectory.py @@ -1,19 +1,22 @@ -import sys import os -import numpy as np -import scipy as sp import shutil +import sys import tempfile +import time as timer import warnings from io import StringIO + +import numpy as np +import pytest +import scipy as sp from scipy import sparse -import time as timer + from mesohops.basis.hops_aux import AuxiliaryVector as AuxVec from mesohops.noise.hops_noise import HopsNoise -from mesohops.util.bath_corr_functions import bcf_convert_sdl_to_exp from mesohops.storage.hops_storage import HopsStorage from mesohops.trajectory.exp_noise import bcf_exp from mesohops.trajectory.hops_trajectory import HopsTrajectory as HOPS +from mesohops.util.bath_corr_functions import bcf_convert_dl_to_exp from mesohops.util.exceptions import UnsupportedRequest from mesohops.util.physical_constants import precision # constant @@ -41,6 +44,18 @@ "PARAM_NOISE1": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]], } +sys_param_noise2 = { + "HAMILTONIAN": np.array([[0, 10.0], [10.0, 0]], dtype=np.float64), + "GW_SYSBATH": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]], + "L_HIER": [loperator[0], loperator[0], loperator[1], loperator[1]], + "L_NOISE1": [loperator[0], loperator[0], loperator[1], loperator[1]], + "L_NOISE2": [loperator[0], loperator[0], loperator[1], loperator[1]], + "ALPHA_NOISE1": bcf_exp, + "ALPHA_NOISE2": bcf_exp, + "PARAM_NOISE1": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]], + "PARAM_NOISE2": [[10.0, 10.0], [5.0, 5.0], [10.0, 10.0], [5.0, 5.0]] +} + hier_param = {"MAXHIER": 4} eom_param = {"TIME_DEPENDENCE": False, "EQUATION_OF_MOTION": "NORMALIZED NONLINEAR"} @@ -155,7 +170,7 @@ def test_initialize(): partial_dict = {key: val for (key, val) in integrator_param.items()} partial_dict['EARLY_INTEGRATOR_STEPS'] = 3 assert hops_partial.integration_param == partial_dict - try: + with pytest.raises(UnsupportedRequest) as excinfo: hops_broken = HOPS( sys_param, noise_param=noise_param, @@ -163,9 +178,8 @@ def test_initialize(): eom_param=eom_param, integration_param=integrator_param_broken, ) - except UnsupportedRequest as excinfo: - assert ("The current code does not support PANCAKE_MIXER in the integration" - in str(excinfo)) + assert ("The current code does not support PANCAKE_MIXER in the integration" + in str(excinfo.value)) # Checks to make sure initialization is timed correctly by measuring control timing # and comparing to the time given by starting a timer, waiting one second, and @@ -399,7 +413,7 @@ def test_inchworm_aux(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([2, 2, 2], dtype=np.float64) gw_sysbath = [] @@ -457,7 +471,7 @@ def test_inchworm_aux(): # First inchworm # ---------------------------------------------------------------------------------- - state_update, aux_update, phi = hops_inchworm.inchworm_integrate( + state_update, aux_update, phi, z_mem = hops_inchworm.inchworm_integrate( state_update, aux_update, 2.0 ) aux_new = aux_update @@ -481,7 +495,7 @@ def test_inchworm_aux(): # Second inchworm # ---------------------------------------------------------------------------------- - state_update, aux_update, phi = hops_inchworm.inchworm_integrate( + state_update, aux_update, phi, z_mem = hops_inchworm.inchworm_integrate( state_update, aux_update, 2.0 ) aux_new = aux_update @@ -528,7 +542,7 @@ def test_inchworm_aux(): # Third inchworm # ---------------------------------------------------------------------------------- - state_update, aux_update, phi = hops_inchworm.inchworm_integrate( + state_update, aux_update, phi, z_mem = hops_inchworm.inchworm_integrate( state_update, aux_update, 2.0 ) aux_new = aux_update @@ -618,7 +632,7 @@ def test_inchworm_state(): e_lambda = 20.0 gamma = 50.0 temp = 140.0 - (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp) + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) loperator = np.zeros([10, 10, 10], dtype=np.float64) gw_sysbath = [] @@ -694,7 +708,7 @@ def test_inchworm_state(): # First inchworm step # ---------------------------------------------------------------------------------- - state_update, aux_update, phi = hops_inchworm.inchworm_integrate( + state_update, aux_update, phi, z_mem = hops_inchworm.inchworm_integrate( state_update, aux_update, 2.0 ) state_new = state_update @@ -703,7 +717,7 @@ def test_inchworm_state(): # Second inchworm step # ---------------------------------------------------------------------------------- - state_update, aux_update, phi = hops_inchworm.inchworm_integrate( + state_update, aux_update, phi, z_mem = hops_inchworm.inchworm_integrate( state_update, aux_update, 2.0 ) state_new = state_update @@ -712,7 +726,7 @@ def test_inchworm_state(): # Third inchworm step # ---------------------------------------------------------------------------------- - state_update, aux_update, phi = hops_inchworm.inchworm_integrate( + state_update, aux_update, phi, z_mem = hops_inchworm.inchworm_integrate( state_update, aux_update, 2.0 ) state_new = state_update @@ -720,6 +734,181 @@ def test_inchworm_state(): assert np.array_equal(state_new, known) +def test_inchworm_z_mem(): + """ + test that noise memory term is being consistently updated during inchworming. + """ + # System parameters and dictionaries + # ---------------------------------- + nsite = 4 + e_lambda = 500.0 + gamma = 50.0 + temp = 300 + (g_0, w_0) = bcf_convert_dl_to_exp(e_lambda, gamma, temp) + + loperator = np.zeros([4, 4, 4], dtype=np.float64) + gw_sysbath = [] + lop_list = [] + for i in range(nsite): + loperator[i, i, i] = 1.0 + gw_sysbath.append([g_0, w_0]) + lop_list.append(sp.sparse.coo_matrix(loperator[i])) + gw_sysbath.append([-1j * np.imag(g_0), 500.0]) + lop_list.append(loperator[i]) + + hs = np.zeros([nsite, nsite], dtype=np.complex128) + hs += np.diag([1000] * (nsite - 1), 1) + np.diag([1000] * (nsite - 1), -1) + + hier_param = {"MAXHIER": 1} + + sys_param = { + "HAMILTONIAN": np.array(hs, dtype=np.complex128), + "GW_SYSBATH": gw_sysbath, + "L_HIER": lop_list, + "L_NOISE1": lop_list, + "ALPHA_NOISE1": bcf_exp, + "PARAM_NOISE1": gw_sysbath, + } + + eom_param = {"EQUATION_OF_MOTION": "NORMALIZED NONLINEAR"} + + # To test the change in z_mem during inchworming, we can use the inchworm cap in + # our parameter dictionaries to limit inchworming and compare what z_mem arrays are + # being produced. The values of z_mem should update with inchworming for the first 2 + # iterations before stabilizing due to the limits of RK4 integration. + # ---------------------------------------------------------------------------------- + # No inchworming + integrator_param_no_inchworm = { + "INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 0, + 'STATIC_BASIS': None + } + + # 1 inchworm iteration + integrator_param_inchworm_1 = { + "INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 1, + 'INCHWORM_CAP': 1, + 'STATIC_BASIS': None + } + + # 2 inchworm iterations + integrator_param_inchworm_2 = { + "INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 1, + 'INCHWORM_CAP': 2, + 'STATIC_BASIS': None + } + + # 3 inchworm iterations + integrator_param_inchworm_3 = { + "INTEGRATOR": "RUNGE_KUTTA", + 'EARLY_ADAPTIVE_INTEGRATOR': 'INCH_WORM', + 'EARLY_INTEGRATOR_STEPS': 1, + 'INCHWORM_CAP': 3, + 'STATIC_BASIS': None + } + + # Initialize and propagate all trajectories by 1 step + # --------------------------------------------------- + psi_0 = np.array([0.0] * nsite, dtype=np.complex128) + psi_0[0] = 1.0 + psi_0 = psi_0 / np.linalg.norm(psi_0) + + hops_no_inchworm = HOPS( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_no_inchworm, + ) + hops_no_inchworm.make_adaptive(1e-15, 1e-15) + hops_no_inchworm.initialize(psi_0) + hops_no_inchworm.propagate(2.0, 2.0) # Performing a single time-step + + hops_inchworm_1 = HOPS( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_inchworm_1, + ) + hops_inchworm_1.make_adaptive(1e-15, 1e-15) + hops_inchworm_1.initialize(psi_0) + hops_inchworm_1.propagate(2.0, 2.0) # Performing a single time-step + + hops_inchworm_2 = HOPS( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_inchworm_2, + ) + hops_inchworm_2.make_adaptive(1e-15, 1e-15) + hops_inchworm_2.initialize(psi_0) + hops_inchworm_2.propagate(2.0, 2.0) # Performing a single time-step + + hops_inchworm_3 = HOPS( + sys_param, + noise_param=noise_param, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_inchworm_3, + ) + hops_inchworm_3.make_adaptive(1e-15, 1e-15) + hops_inchworm_3.initialize(psi_0) + hops_inchworm_3.propagate(2.0, 2.0) # Performing a single time-step + + # When propagating the HOPS wavefunction using the 4th order Runge-Kutta (RK4) + # integrator, the total derivative is calculated over 4 iterations. In each + # iteration, fluxes to adjacent elements in the HOPS wavefunction are permitted. + # Therefore, any modes of a distance less than 4 should be affected by inchworm + # integration. This is represented diagrammatically below, with included basis + # elements listed as Ø, elements outside of the adaptive basis listed as O, and the + # distances determined by counting the linking | or — symbols. As shown below, all + # physical wavefunction elements will be included by the 2nd inchworm step. + # ---------------------------------------------------------------------------------- + # O O O O first auxiliaries for each mode (attached to its respective site) + # | | | | + # O — O — O — O physical wavefunction + # 1-->2-->3-->4 + # + # 0 timesteps + # O O O O + # | | | | + # Ø — O — O — O + # + # 1 timestep, 0 inchworm iterations + # Ø O O O + # | | | | + # Ø — Ø — O — O + # + # 1 timestep, 1 inchworm iteration + # Ø Ø O O + # | | | | + # Ø — Ø — Ø — O + # + # 1 timestep, 2 inchworm iterations + # Ø Ø Ø O + # | | | | + # Ø — Ø — Ø — Ø + # ---------------------------------------------------------------------------------- + + # Testing that z_mem changes on all inchworm iterations that change the physical + # wavefunction basis. There should be at least one significant difference in z_mem + # for each inchworm iteration that adds RK4-accessible physical wavefunction + # elements, but once all physical wavefunction elements are included, inchworming + # should make no further changes to z_mem. + # --------------------------------------------------------------------------------- + # First and second inchworm iterations should change z_mem + assert abs(hops_no_inchworm.z_mem - hops_inchworm_1.z_mem).max() > precision + assert abs(hops_inchworm_1.z_mem - hops_inchworm_2.z_mem).max() > precision + # Final inchworm iteration should not change z_mem + assert abs(hops_inchworm_2.z_mem - hops_inchworm_3.z_mem).max() <= precision + def test_prepare_zstep(): """ @@ -802,12 +991,10 @@ def test_early_time_integrator(): ) hops.make_adaptive(0.001, 0.001) hops.initialize(psi_0) - try: + with pytest.raises(UnsupportedRequest) as excinfo: hops.propagate(2.0, 2.0) - except UnsupportedRequest as excinfo: - if "does not support CHOOSE_BASIS_RANDOMLY in the early time integrator " \ - "clause" not in str(excinfo): - pytest.fail() + assert "does not support CHOOSE_BASIS_RANDOMLY in the early time integrator " \ + "clause" in str(excinfo.value) def test_static_state_inchworm_hierarchy(): """ @@ -1456,3 +1643,247 @@ def test_save_slices_mixed_data_types(): # Clean up shutil.rmtree(temp_dir) +def test_initialize_2_noise(capsys): + """ + Tests that riskily initializing the second noise prints the correct warnings. + """ + noise2_param_seed_clash = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + } + + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_seed_clash, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert any("Using the same seed for both noise 1 and" in str(warning.message) + for warning in w) + + noise2_param_correct = { + "SEED": 1010101, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + } + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_correct, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not any("Using the same seed for both noise 1 and" in str( + warning.message) for warning in w) + + noise2_param_none = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + } + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_none, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not any("Using the same seed for both noise 1 and" in str( + warning.message) for warning in w) + # Tests that when the time axes are identical, we don't warn for mismatched + # time axes. + assert not any("Time axes of noise 1 and noise 2 are" in str( + warning.message) for warning in w) + + noise_param_none = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + } + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param_none, + noise2_param=noise2_param_none, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not any("Using the same seed for both noise 1 and" in str( + warning.message) for warning in w) + + + noise2_param_mismatch_tlen = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 20.0, # Units: fs + "TAU": 1.0, # Units: fs + } + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_mismatch_tlen, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert any("Time axes of noise 1 and noise 2 are" in str( + warning.message) for warning in w) + + # Also tests explicitly false interpolation + noise2_param_mismatch_tau = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 0.5, # Units: fs + "INTERPOLATE": False + } + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_mismatch_tau, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert any("Time axes of noise 1 and noise 2 are" in str( + warning.message) for warning in w) + + noise2_param_interp = { + "SEED": None, + "MODEL": "FFT_FILTER", + "TLEN": 20.0, # Units: fs + "TAU": 0.5, # Units: fs + "INTERPOLATE": True + } + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_interp, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not any("Time axes of noise 1 and noise 2 are" in str( + warning.message) for warning in w) + + # Check that HopsTrajectory and subsidiary helper functions manage FLAG_REAL + # correctly. These are technically integrated tests because the result should be + # agnostic to the handling: noise1 should force FLAG_REAL to be False and noise2 + # should assume FLAG_REAL is True, but allow it to be manually set False. + noise_param_real = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + "FLAG_REAL": True, + } + + noise2_param_real = { + "SEED": 1010101, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + "FLAG_REAL": True, + } + + noise_param_complex = { + "SEED": 0, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + "FLAG_REAL": False, + } + + noise2_param_complex = { + "SEED": 1010101, + "MODEL": "FFT_FILTER", + "TLEN": 10.0, # Units: fs + "TAU": 1.0, # Units: fs + "FLAG_REAL": False, + } + + # If noise 1 is flagged real, then a warning will be raised. Similarly, + # if FLAG_REAL for noise 2 is not specified, a warning will be raised to let the + # user know that it has been defaulted True. + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param_real, + noise2_param=noise2_param_real, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert any("Noise 1 should never be flagged real" in str( + warning.message) for warning in w) + assert not any("Noise 2 FLAG_REAL not specified: setting to True." in str( + warning.message) for warning in w) + + with warnings.catch_warnings(record=True) as w: + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_correct, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not any("Noise 1 should never be flagged real" in str( + warning.message) for warning in w) + assert any("Noise 2 FLAG_REAL not specified: setting to True." in str( + warning.message) for warning in w) + + + # HopsTrajectory should force noise1 to not be FLAG_REAL + hops = HOPS( + sys_param_noise2, + noise_param=noise_param_real, + noise2_param=noise2_param_real, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not hops.noise1.param["FLAG_REAL"] + assert hops.noise2.param["FLAG_REAL"] + + # HopsTrajectory should default noise1 to not be FLAG_REAL and noise2 to be + # FLAG_REAL + hops = HOPS( + sys_param_noise2, + noise_param=noise_param, + noise2_param=noise2_param_correct, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not hops.noise1.param["FLAG_REAL"] + assert hops.noise2.param["FLAG_REAL"] + + # Explicitly setting FLAG_REAL to False should always work for both noise1 and + # noise2. + hops = HOPS( + sys_param_noise2, + noise_param=noise_param_complex, + noise2_param=noise2_param_complex, + hierarchy_param=hier_param, + eom_param=eom_param, + integration_param=integrator_param_empty, + ) + assert not hops.noise1.param["FLAG_REAL"] + assert not hops.noise2.param["FLAG_REAL"] \ No newline at end of file diff --git a/tests/test_noise_fft.py b/tests/test_noise_fft.py index 2d0d42b..5690349 100644 --- a/tests/test_noise_fft.py +++ b/tests/test_noise_fft.py @@ -1,10 +1,12 @@ import os + import numpy as np -from mesohops.trajectory.exp_noise import bcf_exp +import pytest + from mesohops.noise.hops_noise import HopsNoise +from mesohops.trajectory.exp_noise import bcf_exp from mesohops.util.exceptions import UnsupportedRequest - __title__ = "Test of FFT_FILTER noise model" __author__ = "J. K. Lynd" __version__ = "1.2" @@ -357,24 +359,19 @@ def test_prepare_rand(): noise_param["SEED"] = np.array([np.arange(17), np.arange(17)]) test_noise = HopsNoise(noise_param, noise_corr) - try: - test_rand_improper_list_seed = test_noise._prepare_rand() - except UnsupportedRequest as excinfo: - assert "Noise.param[SEED] is an array of the wrong length" in str(excinfo) + with pytest.raises(UnsupportedRequest, match="Noise.param\[SEED\] is an array of " + "the wrong length"): + test_noise._prepare_rand() noise_param["SEED"] = "string" test_noise = HopsNoise(noise_param, noise_corr) - try: - test_rand_string_seed = test_noise._prepare_rand() - except UnsupportedRequest as excinfo: - assert "is not the address of a valid file" in str(excinfo) + with pytest.raises(UnsupportedRequest, match="is not the address of a valid file"): + test_noise._prepare_rand() - try: + with pytest.raises(TypeError, match='is of type'): noise_param['SEED'] = HopsNoise({}, {"N_L2": 1}) test_noise = HopsNoise(noise_param, noise_corr) - test_rand_HopsNoise = test_noise._prepare_rand() - except TypeError as excinfo: - assert 'is of type' in str(excinfo) + test_noise._prepare_rand() def test_construct_indexing(): """ diff --git a/tests/test_noise_zero.py b/tests/test_noise_zero.py index 6b1958f..34056f7 100644 --- a/tests/test_noise_zero.py +++ b/tests/test_noise_zero.py @@ -88,4 +88,7 @@ def test_noiseModel(): assert np.allclose(testnoise1, np.zeros_like(testnoise1)) assert np.size(noiseModel.get_noise(noiseModel.param["T_AXIS"])) == np.size( noiseModel.param["T_AXIS"])*noiseModel.param["N_L2"] + # The stored noise should be an integer 0 - it is never actually accessed, + # so a simple-to-understand placeholder value is used instead. + assert noiseModel._noise == 0 diff --git a/tests/utils.py b/tests/utils.py new file mode 100644 index 0000000..bd4729d --- /dev/null +++ b/tests/utils.py @@ -0,0 +1,81 @@ +# utils.py +# Utility functions for testing in the mesohops library. +# Provides deep comparison utilities for complex data structures such as numpy arrays, sparse matrices, lists, and dictionaries. + +import numpy as np # Import numpy for array operations and testing utilities +import scipy as sp # Import scipy for sparse matrix support + + +def _compare_values(v1, v2): + """ + Recursively compares two values for equality, supporting a variety of data types. + Handles numpy arrays, scipy sparse matrices, lists, dictionaries, and scalars. + Raises an AssertionError if the values are not considered equal. + + Inputs + ------ + 1. v1 : Any + First value to compare. + 2. v2 : Any + Second value to compare. + + Returns + ------- + None + """ + # Compare sparse matrices by converting to dense arrays + if sp.sparse.issparse(v1) or sp.sparse.issparse(v2): + # Both must be sparse matrices of compatible shape and type + assert sp.sparse.issparse(v1) and sp.sparse.issparse(v2), "Both values must be sparse matrices." + np.testing.assert_allclose(v1.toarray(), v2.toarray()) + + # Compare numpy arrays + elif isinstance(v1, np.ndarray) or isinstance(v2, np.ndarray): + # Handle object arrays element-wise (recursively) + if isinstance(v1, np.ndarray) and v1.dtype == object: + assert isinstance(v2, np.ndarray) and v2.dtype == object, "Both arrays must be object dtype." + assert v1.shape == v2.shape, "Object arrays must have the same shape." + for a, b in zip(v1.flat, v2.flat): + _compare_values(a, b) + else: + # Use numpy's allclose for numerical arrays + np.testing.assert_allclose(v1, v2) + + # Recursively compare lists + elif isinstance(v1, list) and isinstance(v2, list): + assert len(v1) == len(v2), "Lists must have the same length." + for a, b in zip(v1, v2): + _compare_values(a, b) + + # Recursively compare dictionaries + elif isinstance(v1, dict) and isinstance(v2, dict): + assert set(v1.keys()) == set(v2.keys()), "Dictionaries must have the same keys." + for k in v1: + _compare_values(v1[k], v2[k]) + + else: + # Fallback for scalars and other comparable objects + assert v1 == v2, f"Values {v1} and {v2} are not equal." + + +def compare_dictionaries(dict1, dict2): + """ + Assert that two parameter dictionaries are equivalent, deeply comparing all values. + Useful for testing that two sets of parameters or results are identical in structure and content. + + Parameters + ---------- + 1. dict1 : dict + First dictionary for comparison. + 2. dict2 : dict + Second dictionary for comparison. + + Returns + ------- + None + """ + # Ensure both dictionaries have the same set of keys + assert set(dict1.keys()) == set(dict2.keys()), "Dictionaries must have the same keys." + for key in dict1: + print(f'Comparing {key}') # Print the key being compared for debugging + _compare_values(dict1[key], dict2[key])