import sys import traceback import copy import ctypes as c import numpy as np import os import mpmath def _print_error(msg_str): sys.stderr.write(msg_str + "\n") sys.stderr.flush() try: _OpenBLASLibPath = os.path.join(os.path.dirname(os.path.abspath(__file__)),"OpenBLAS_install/lib/libopenblas.so.0") c.CDLL(_OpenBLASLibPath) except OSError: _print_error("Warning: Failed to load OpenBLAS in its default place. " "Excepting the library to be available from the system.") _SpintrumLibPath = os.path.join(os.path.dirname(os.path.abspath(__file__)), "lib/libspintrum.so") if not os.path.exists(_SpintrumLibPath): msg = "WARNING: Could not load the spintrum shared library in path " \ + _SpintrumLibPath + ". This means that the installation is not valid " \ "and will most likely not work." _print_error(msg) print(msg) mpmath.mp.dps = 50 def fft_spectral_density(data_list, sampling_rate): """ Calculates the amplitude spectral density with both x and y axes :param data_list: Data, for which the FFT to be taken :param sampling_rate: the sample rate of the data :return: a dict with two elements, "x" for the x-axis (frequency axis), and "y" for the y-axis (the amplitude axis with units Amplitude/sqrt(Hz)) """ N = len(data_list) fft_data = np.fft.rfft(data_list)[0:int(N / 2 + 1)] bandwidth = np.double(sampling_rate / 2) scaling_factor = np.sqrt(1/bandwidth/np.double(N)) x_axis = sampling_rate * np.array(range(int(N / 2 + 1))) / (np.double(N)) abs_fft_data = scaling_factor*np.array(list(map(np.absolute,fft_data))) return {"x": x_axis, "y": abs_fft_data} def FFTSpectralDensity(data_list, sampling_rate): """ Deprecated. Please use fft_spectral_density() """ return fft_spectral_density(data_list, sampling_rate) def who_am_i(): stack = traceback.extract_stack() filename, code_line, func_name, text = stack[-2] return func_name def _raise_error(msg_str): sys.stderr.write(msg_str + "\n") raise Exception(msg_str) def _is_var_type_correct(value, type_or_types): if type(type_or_types) is list: for tp in type_or_types: if type(value) is tp: return True return False else: if type(value) is type_or_types: return True else: return False def _get_types_list_types_names(type_or_types): if type(type_or_types) is list: return str([tp.__name__ for tp in type_or_types]) else: return str(type_or_types.__name__) def _check_dict_elements_types(input_dict, types_dict): for val in input_dict: if val in types_dict: if _is_var_type_correct(input_dict[val], types_dict[val]) is False: _raise_error("Error: Type is not of expected type for variable: " + str(val) + ". Possible allowed type(s) are: " + _get_types_list_types_names(types_dict[val])) def _check_elements_exist_in_dict(input_dict, must_exist, can_exist): for el in must_exist: if el not in input_dict: _raise_error("Error: You should specify a value for: " + el) joined_elements = must_exist + can_exist for el in input_dict: if el not in joined_elements: _raise_error("Error: You have specified an unknown variable: " + el) def _analyze_direction_string(dirs): activated_directions = {'x': 0, 'y': 0, 'z': 0, 'dirs': dirs} if (len(dirs) != 0) and ("x" not in dirs) and ("y" not in dirs) and ("z" not in dirs): _raise_error("Directions, for which trace has to be stored must be given as a string; eg x or xyz or xz, etc.") else: for c in dirs: if c == 'x': activated_directions['x'] = 1 elif c == 'y': activated_directions['y'] = 1 elif c == 'z': activated_directions['z'] = 1 return activated_directions class SpinOperations: _OPERATION_NAME_IDENTIFIER = '_OperationIdentifier' OPERATION__THERMAL_POPULATE = 0 OPERATION__TIP_SPINS = 1 OPERATION__SET_HAMILTONIAN = 2 OPERATION__EVOLVE_TIME_INDEPENDENT = 3 OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION = 4 OPERATION__EVOLVE_TIME_DEPENDENT = 5 def __init__(self): self.operations_stack = None self.total_num_of_points = 0 self.clear_operations() def clear_operations(self): self.operations_stack = [] def add_operation(self, operation, parameters={}): if type(parameters) is not dict: _raise_error(who_am_i() + ": Error: parameters should be a dict.") if operation == SpinOperations.OPERATION__THERMAL_POPULATE: argument_types = {'Bx': [int, float], 'By': [int, float], 'Bz': [int, float],'T': [int,float]} must_exist_arguments = ['Bx', 'By', 'Bz','T'] optional_arguments = [] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) params_copy = copy.deepcopy(parameters) params_copy.update({SpinOperations._OPERATION_NAME_IDENTIFIER: operation}) self.operations_stack.append(params_copy) elif operation == SpinOperations.OPERATION__TIP_SPINS: argument_types = {'direction': str, 'BVsTArea': [int, float]} must_exist_arguments = ['direction', 'BVsTArea'] optional_arguments = [] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) params_copy = copy.deepcopy(parameters) params_copy.update({SpinOperations._OPERATION_NAME_IDENTIFIER: operation}) self.operations_stack.append(params_copy) elif operation == SpinOperations.OPERATION__SET_HAMILTONIAN: argument_types = {'Bx': [int, float], 'By': [int, float], 'Bz': [int, float], 'samplingRate': [int, float]} must_exist_arguments = [] optional_arguments = ['Bx', 'By', 'Bz'] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) params_copy = copy.deepcopy(parameters) params_copy.update({SpinOperations._OPERATION_NAME_IDENTIFIER: operation}) self.operations_stack.append(params_copy) elif operation == SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT: argument_types = {'points': int, 'threads': int} must_exist_arguments = ['points', 'threads'] optional_arguments = [] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) params_copy = copy.deepcopy(parameters) params_copy.update({SpinOperations._OPERATION_NAME_IDENTIFIER: operation}) self.total_num_of_points += parameters['points'] self.operations_stack.append(params_copy) elif operation == SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION: argument_types = {'samplingRate': [int, float], 'measurementDirection': str} must_exist_arguments = ['samplingRate', 'measurementDirection'] optional_arguments = [] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) params_copy = copy.deepcopy(parameters) params_copy.update({SpinOperations._OPERATION_NAME_IDENTIFIER: operation}) self.operations_stack.append(params_copy) elif operation == SpinOperations.OPERATION__EVOLVE_TIME_DEPENDENT: argument_types = {'threads': int, 'input': dict, 'samplingRate': [int, float], 'store': str} must_exist_arguments = ['threads', 'input', 'samplingRate', ] optional_arguments = ['store'] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) # use set to remove duplicates, and restore in dirs_to_store to get the reduced version # (removing duplicates preserves order by using key=parameters['store'].index) parameters['dirs_to_store'] = \ _analyze_direction_string(''.join(sorted(set(parameters['store']), key=parameters['store'].index))) # total number of points = number of input values*number of directions to measure self.total_num_of_points += \ len(list(parameters['input'].values())[0])*len(parameters['dirs_to_store']['dirs']) params_copy = copy.deepcopy(parameters) params_copy.update({SpinOperations._OPERATION_NAME_IDENTIFIER: operation}) self.operations_stack.append(params_copy) else: _raise_error("Error: Operation Unknown!") def _check_j_couplings_sanity(j_couplings): """ raises an error if j couplings don't make sense :param j_couplings: :return: """ if j_couplings is None: #no j-couplings provided return for i in range(len(j_couplings)): if len(j_couplings[i]) != len(j_couplings): _raise_error("Error: j-couplings must be a square matrix.") for j in range(len(j_couplings[i])): if j_couplings[i][j] != 0 and i >= j: _raise_error("Error: j-couplings must be an upper triangular matrix. " "All other elements must be exactly 0.") def simulate(**kwargs): argument_types = {'gyromagneticRatios': list, 'jCouplings': list, 'spinMultiplicities': list, 'spinOperations': SpinOperations} must_exist_arguments = ['gyromagneticRatios', 'spinOperations'] optional_arguments = ['jCouplings','spinMultiplicities'] _check_elements_exist_in_dict(kwargs, must_exist_arguments, optional_arguments) _check_dict_elements_types(kwargs, argument_types) if 'jCouplings' not in kwargs: kwargs['jCouplings'] = None if 'spinMultiplicities' not in kwargs: kwargs['spinMultiplicities'] = None _check_j_couplings_sanity(kwargs['jCouplings']) if kwargs['spinMultiplicities'] is not None and \ len(kwargs['gyromagneticRatios']) != len(kwargs['spinMultiplicities']): _raise_error("Spin multiplicities must have a length equal to the gyromagnetic ratios, " "or alternatively, if you keep it empty, then it is assumed that all spins are 1/2") # print(kwargs) return _get_signal_c_call(kwargs) class _CSpinOperation(c.Structure): _fields_ = [("OperationIdentifier", c.c_long), ("Params_int_count", c.c_long), ("Params_int", c.POINTER(c.c_long)), ("Params_real_count", c.c_long), ("Params_real", c.POINTER(c.c_double)), ("Params_array_count", c.c_long), ("Params_array_lengths", c.POINTER(c.c_long)), ("Params_array_names", c.POINTER(c.POINTER(c.c_char))), ("Params_array", c.POINTER(c.POINTER(c.c_double))) ] def _direction_str_to_int(dir_str): if dir_str.lower() == 'x': return 0 elif dir_str.lower() == 'y': return 1 elif dir_str.lower() == 'z': return 2 else: _raise_error("Direction of the tipping magnetic field can be either x, y or z. " "You've entered \"" + dir_str + "\"") def _create_c_operation(operation_id, int_params, real_params, array_params): # if int params exist, make them in a pointer to an array, otherwise, make it NULL (None) if len(int_params) > 0: int_params_c = (c.c_long * len(int_params))(*int_params) #convert list to c-array int_params_ptr = c.cast(int_params_c, c.POINTER(c.c_long)) #cast array to pointer else: int_params_ptr = None #if real params exist, make them in a pointer to an array, otherwise, make it NULL (None) if len(real_params) > 0: real_params_c = (c.c_double * len(real_params))(*real_params) #convert list to c-array real_params_ptr = c.cast(real_params_c, c.POINTER(c.c_double)) #cast array to pointer else: real_params_ptr = None if len(array_params) > 0: # arrays lengths array_lengths_list = list(map(len,array_params.values())) data_array_lengths_c = (c.c_long * len(array_params))(*array_lengths_list) #convert list to c-array data_array_lengths_ptr = c.cast(data_array_lengths_c, c.POINTER(c.c_long)) #cast array to pointer # arrays labels labels_array_type = (c.POINTER(c.c_char)*len(array_params)) labels_c_array = labels_array_type() labels_keys_list = list(array_params.keys()) for i in range(len(labels_keys_list)): labels_c_array[i] = c.create_string_buffer(str.encode(labels_keys_list[i])) # arrays data arrays arrays_array_type = (c.POINTER(c.c_double) * len(array_params)) arrays_c_array = arrays_array_type() arrays_list = list(array_params.values()) for i in range(len(arrays_list)): single_array_c_type = (c.c_double * len(arrays_list[i]))(*arrays_list[i]) # convert list to c-array single_array_c = single_array_c_type arrays_c_array[i] = c.cast(single_array_c, c.POINTER(c.c_double)) #cast array to pointer else: data_array_lengths_ptr = None labels_c_array = None arrays_c_array = None op_values = _CSpinOperation(operation_id, len(int_params), int_params_ptr, len(real_params), real_params_ptr, len(array_params), data_array_lengths_ptr, labels_c_array, arrays_c_array) return c.byref(op_values) def _get_operations_pointers(operations_stack): num_of_ops = len(operations_stack) structs_ptr_array_type = c.c_void_p * num_of_ops structs_ptr_array_obj = structs_ptr_array_type() structs_array = [] for op in operations_stack: if op["_OperationIdentifier"] == SpinOperations.OPERATION__THERMAL_POPULATE: structs_array.append(_create_c_operation(op["_OperationIdentifier"], [], [op['Bx'],op['By'],op['Bz'],op['T']], {})) elif op["_OperationIdentifier"] == SpinOperations.OPERATION__TIP_SPINS: structs_array.append(_create_c_operation(op["_OperationIdentifier"],[_direction_str_to_int(op['direction'])], [op['BVsTArea']], {})) elif op["_OperationIdentifier"] == SpinOperations.OPERATION__SET_HAMILTONIAN: structs_array.append(_create_c_operation(op["_OperationIdentifier"], [], [op["Bx"], op["By"], op["Bz"]], {})) elif op["_OperationIdentifier"] == SpinOperations.OPERATION__EVOLVE_TIME_INDEPENDENT: structs_array.append(_create_c_operation(op["_OperationIdentifier"], [op["points"],op["threads"]], [], {})) elif op["_OperationIdentifier"] == SpinOperations.OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION: structs_array.append(_create_c_operation(op["_OperationIdentifier"], [_direction_str_to_int(op['measurementDirection'])], [op["samplingRate"]], {})) elif op["_OperationIdentifier"] == SpinOperations.OPERATION__EVOLVE_TIME_DEPENDENT: structs_array.append(_create_c_operation(op["_OperationIdentifier"], [op["threads"], op['dirs_to_store']['x'], op['dirs_to_store']['y'], op['dirs_to_store']['z']], [op['samplingRate']], op["input"])) return structs_array def _get_signal_c_call(parameters): structs_array = _get_operations_pointers(parameters['spinOperations'].operations_stack) operations_stack = parameters['spinOperations'].operations_stack num_of_ops = len(operations_stack) structs_ptr_array_type = c.c_void_p * num_of_ops structs_ptr_array_obj = structs_ptr_array_type() # print(structs_ptr_array_obj) # This next loop cannot be included in the function _get_operations_pointers # because taking c-pointers to structs_array elements doesn't make python preserve them, which causes them to # be deleted when exiting the function for i in range(len(structs_array)): structs_ptr_array_obj[i] = c.cast(structs_array[i], c.c_void_p) gammas = parameters['gyromagneticRatios'] mults = parameters['spinMultiplicities'] j_couplings = _convert_j_couplings_to_1d_array(parameters['jCouplings']) gammas_c_version = (c.c_double * len(gammas))(*gammas) gammas_ptr = c.cast(gammas_c_version, c.POINTER(c.c_double * len(gammas))) if mults is None: mults_ptr = c.c_void_p(None) else: mults_c_version = (c.c_int * len(mults))(*mults) mults_ptr = c.cast(mults_c_version, c.POINTER(c.c_int * len(mults))) if j_couplings == []: j_couplings_ptr = c.c_void_p(None) else: j_couplings_c_version = (c.c_double * len(j_couplings))(*j_couplings) j_couplings_ptr = c.cast(j_couplings_c_version, c.POINTER(c.c_double * len(j_couplings))) mylib = c.CDLL(_SpintrumLibPath) num_of_points = parameters['spinOperations'].total_num_of_points return_array_type = c.c_double * num_of_points returned_array_ptr_type = c.POINTER(return_array_type) return_array = return_array_type() returned_array_ptr = c.byref(return_array) prototype = c.CFUNCTYPE( c.c_long, # return type c.c_void_p, c.c_void_p, c.c_void_p, c.c_void_p, c.c_int, c.c_int, c.c_long, returned_array_ptr_type ) create_signal_function = prototype(('SimulateSignal', mylib)) res = create_signal_function(structs_ptr_array_obj, gammas_ptr, j_couplings_ptr, mults_ptr, c.c_int(len(gammas)), c.c_int(len(structs_ptr_array_obj)), c.c_long(num_of_points), returned_array_ptr) # print(res) output = np.array(c.cast(returned_array_ptr,returned_array_ptr_type).contents) # print(output) return output def _convert_j_couplings_to_1d_array(j_couplings): if j_couplings is None: return [] output_array = [] for i in range(len(j_couplings)): output_array += j_couplings[i][i + 1:len(j_couplings)] return output_array def _get_derivative(params, param_index, fitted_function, epsilon=np.sqrt(np.finfo(float).eps)): """ Calculates the derivatives of the fitted function (after having fitted a function) to some model :param params: list of parameter values at the point where the derivative has to be calculated :param param_index: index of the parameter, for which the derivative has to be taken :param fitted_function: the function, for which the derivative has to be calculated, where this function has to take "params" as an argument :param epsilon: the small step to calculate the numerical difference; the default value uses Python floats' epislon :return: the derivative value of the i-th parameter in params in the function fitted_function """ h = params[param_index] * epsilon low_point__1 = np.delete(np.insert(params, param_index, params[param_index] - h), param_index + 1) high_point_1 = np.delete(np.insert(params, param_index, params[param_index] + h), param_index + 1) low_point__2 = np.delete(np.insert(params, param_index, params[param_index] - 2 * h), param_index + 1) high_point_2 = np.delete(np.insert(params, param_index, params[param_index] + 2 * h), param_index + 1) low_point__3 = np.delete(np.insert(params, param_index, params[param_index] - 3 * h), param_index + 1) high_point_3 = np.delete(np.insert(params, param_index, params[param_index] + 3 * h), param_index + 1) func_value_lp1 = fitted_function(low_point__1) func_value_hp1 = fitted_function(high_point_1) func_value_lp2 = fitted_function(low_point__2) func_value_hp2 = fitted_function(high_point_2) func_value_lp3 = fitted_function(low_point__3) func_value_hp3 = fitted_function(high_point_3) derivative = 1.5*(func_value_hp1['y'] - func_value_lp1['y'])/(2*h) - \ 0.6*(func_value_hp2['y'] - func_value_lp2['y'])/(4*h) + \ 0.1*(func_value_hp3['y'] - func_value_lp3['y'])/(6*h) return derivative def get_correlation_matrix(optimum_params, fitted_function, objective_function): """ A function that calculates the correlation matrix according to the book: Wolberg 2005: Data analysis using least-squares method :param optimum_params: list of optimum fitted parameters :param fitted_function: function that is fitted, that takes optimum_params as a parameter :param objective_function: the objective function used to fit. It should take optimum_params as parameter :return: the correlation matrix calculated """ n = len(optimum_params) # number of parameters p = len(fitted_function(optimum_params)['x']) # number of fitted points min_obj_func_val = objective_function(optimum_params) # minimum objective_function value derivatives_matrix = \ [_get_derivative(optimum_params, i, fitted_function).tolist() for i in range(len(optimum_params))] derivatives_matrix_high_precision = mpmath.matrix(derivatives_matrix) c_matrix = derivatives_matrix_high_precision*derivatives_matrix_high_precision.T factor = (min_obj_func_val**2/(p-n)) correlation_matrix = (c_matrix**-1)*factor return correlation_matrix def get_errorbars(matrix): error_bars = np.array([]) for i in range(matrix.cols): error_bars = np.append(error_bars, np.float(np.sqrt(matrix[i,i]))) return error_bars class AtomicSpecies: GyromagneticRatio = None Name = None @staticmethod def _get_species_class_list(): return AtomicSpecies.__subclasses__() @staticmethod def get_species_names(): return list([cls.Name for cls in AtomicSpecies._get_species_class_list()]) @staticmethod def get_species_gyromagnetic_ratios(): return dict({cls.Name: cls.GyromagneticRatio for cls in AtomicSpecies._get_species_class_list()}) @staticmethod def get_species_list_gyromagnetic_ratios(listOfNames): gyromagnetic_ratios = dict({cls.Name: cls.GyromagneticRatio for cls in AtomicSpecies._get_species_class_list()}) output = [] for atom in listOfNames: output.append(gyromagnetic_ratios[atom]) return output def __init(self): pass class AtomC13(AtomicSpecies): Name = "C13" GyromagneticRatio = 1070.8 def __init__(self): AtomicSpecies.__init__(self) class AtomH1(AtomicSpecies): Name = "H1" GyromagneticRatio = 4257.7 def __init__(self): AtomicSpecies.__init__(self) class SpinSimulator: def __init__(self, **kwargs): self._is_initialized = False self._spin_operations = None self.init(kwargs) def init(self,parameters): # pointer to object self.c_spin_ptr_to_object_type = c.c_void_p # pointer to pointer to object. Necessary to modify the pointer to object self.c_spin_ptr_to_ptr_type = c.POINTER(self.c_spin_ptr_to_object_type) self.c_spin_ptr_to_object = self.c_spin_ptr_to_object_type() self.c_spin_ptr_to_ptr = c.byref(self.c_spin_ptr_to_object) argument_types = {'gyromagneticRatios': list, 'jCouplings': list, 'doPrint': bool} must_exist_arguments = ['gyromagneticRatios'] optional_arguments = ['jCouplings', 'spinMultiplicities', 'doPrint'] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) if 'jCouplings' not in parameters: parameters['jCouplings'] = None if 'spinMultiplicities' not in parameters: parameters['spinMultiplicities'] = None _check_j_couplings_sanity(parameters['jCouplings']) if parameters['spinMultiplicities'] is not None and \ len(parameters['gyromagneticRatios']) != len(parameters['spinMultiplicities']): _raise_error("Spin multiplicities must have a length equal to the gyromagnetic ratios, " "or alternatively, if you keep it empty, then it is assumed that all spins are 1/2") gammas = parameters['gyromagneticRatios'] mults = parameters['spinMultiplicities'] jCouplings = _convert_j_couplings_to_1d_array(parameters['jCouplings']) gammas_c_version = (c.c_double * len(gammas))(*gammas) gammas_ptr = c.cast(gammas_c_version, c.POINTER(c.c_double * len(gammas))) if mults is None: mults_ptr = c.c_void_p(None) else: mults_c_version = (c.c_int * len(mults))(*mults) mults_ptr = c.cast(mults_c_version, c.POINTER(c.c_int * len(mults))) if jCouplings == []: jCouplings_ptr = c.c_void_p(None) else: jCouplings_c_version = (c.c_double * len(jCouplings))(*jCouplings) jCouplings_ptr = c.cast(jCouplings_c_version, c.POINTER(c.c_double * len(jCouplings))) mylib = c.CDLL(_SpintrumLibPath) prototype = c.CFUNCTYPE( c.c_int, # return type c.c_void_p, c.c_void_p, c.c_void_p, c.c_int, c.c_bool, self.c_spin_ptr_to_ptr_type ) initSimulatorCreatorFunction = prototype(('InitializeSimulator', mylib)) res = initSimulatorCreatorFunction(gammas_ptr, jCouplings_ptr, mults_ptr, c.c_int(len(gammas)), parameters['doPrint'], self.c_spin_ptr_to_ptr) if res >= 0: self._is_initialized = True else: self._is_initialized = False _raise_error("An error occurred while initializing the spin simulator. Possibly a memory issue.") return res def update_parameters(self,**parameters): argument_types = {'gyromagneticRatios': list, 'jCouplings': list, 'spinOperations': SpinOperations} must_exist_arguments = [] optional_arguments = ['gyromagneticRatios', 'spinOperations', 'jCouplings'] _check_elements_exist_in_dict(parameters, must_exist_arguments, optional_arguments) _check_dict_elements_types(parameters, argument_types) there_are_gammas = False there_are_j_couplings = False there_are_spin_operations = False if 'gyromagneticRatios' in parameters: gammas = parameters['gyromagneticRatios'] gammas_c_version = (c.c_double * len(gammas))(*gammas) gammas_ptr = c.cast(gammas_c_version, c.POINTER(c.c_double * len(gammas))) there_are_gammas = True if 'jCouplings' in parameters: _check_j_couplings_sanity(parameters['jCouplings']) jCouplings = _convert_j_couplings_to_1d_array(parameters['jCouplings']) jCouplings_c_version = (c.c_double * len(jCouplings))(*jCouplings) jCouplings_ptr = c.cast(jCouplings_c_version, c.POINTER(c.c_double * len(jCouplings))) there_are_j_couplings = True if 'spinOperations' in parameters: spin_ops = parameters['spinOperations'] there_are_spin_operations = True if there_are_gammas or there_are_j_couplings: mylib = c.CDLL(_SpintrumLibPath) prototype = c.CFUNCTYPE( c.c_int, # return type c.c_void_p, c.c_void_p ) #set gammas values if there_are_gammas: set_gammas_function = prototype(('SetGammasInSimulator', mylib)) res = set_gammas_function(gammas_ptr, self.c_spin_ptr_to_object) if res >= 0: pass else: _raise_error("An error was found while trying to set gamma values.") # set j-couplings values if there_are_j_couplings: set_j_couplings_function = prototype(('SetJCouplingsInSimulator', mylib)) res = set_j_couplings_function(jCouplings_ptr, self.c_spin_ptr_to_object) if res >= 0: pass else: _raise_error("An error was found while trying to set j-couplings values.") if there_are_spin_operations: self._spin_operations = spin_ops def simulate(self): if self._spin_operations is None: _raise_error("SpinOperations not initialized. Use update_parameters() and set spinOperations.") structs_array = _get_operations_pointers(self._spin_operations.operations_stack) num_of_ops = len(self._spin_operations.operations_stack) structs_ptr_array_type = c.c_void_p * num_of_ops structs_ptr_array_obj = structs_ptr_array_type() for i in range(len(structs_array)): structs_ptr_array_obj[i] = c.cast(structs_array[i], c.c_void_p) mylib = c.CDLL(_SpintrumLibPath) num_of_points = self._spin_operations.total_num_of_points returnArrayType = c.c_double * num_of_points returnedArrayPtrType = c.POINTER(returnArrayType) returnArray = returnArrayType() returnedArrayPtr = c.byref(returnArray) prototype = c.CFUNCTYPE( c.c_long, #return type self.c_spin_ptr_to_object_type, c.c_void_p, c.c_int, c.c_long, returnedArrayPtrType ) createSpectFunction = prototype(('CreateSignalFromInitializedSimulator', mylib)) res = createSpectFunction(self.c_spin_ptr_to_object, structs_ptr_array_obj, c.c_int(len(structs_ptr_array_obj)), c.c_long(num_of_points), returnedArrayPtr) # print(res) output = np.array(c.cast(returnedArrayPtr,returnedArrayPtrType).contents) # print(output) return output def __del__(self): # print("System initialization state: " + str(self._is_initialized)) if self._is_initialized: mylib = c.CDLL(_SpintrumLibPath) prototype = c.CFUNCTYPE( c.c_int, # return type c.c_void_p ) deinitSimulatorFunction = prototype(('DenitializeSimulator', mylib)) res = deinitSimulatorFunction(self.c_spin_ptr_to_object) if res >= 0: self._is_initialized = False else: _raise_error("An error was found while trying to free the memory.")