Spintrum/spintrum/spintrum.py

748 lines
31 KiB
Python
Raw Normal View History

import sys
import traceback
import copy
import ctypes as c
import numpy as np
import os
import mpmath
_SpintrumLibPath = os.path.join(os.path.dirname(os.path.abspath(__file__)),"lib/libspintrum.so")
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 _print_error(msg_str):
sys.stderr.write(msg_str + "\n")
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)
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.")