Spintrum/spintrum/spintrum.py

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")


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.")