diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index fdd3c6dfc..77547d49b 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -91,7 +91,10 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install macos system dependencies - run: brew install openblas fftw + run: | + brew install openblas fftw + which clang + clang --version - name: Install dependencies run: | @@ -125,7 +128,10 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install macos system dependencies - run: brew install openblas fftw + run: | + brew install openblas fftw + which clang + clang --version - name: Install dependencies run: | diff --git a/setup.py b/setup.py index 80d1be385..d4ff8fc7b 100644 --- a/setup.py +++ b/setup.py @@ -66,7 +66,7 @@ def check_if_header_exists(self, header): def conda_setup_for_windows(self): self.libraries += ["fftw3", "openblas"] - self.extra_compile_args = ["/DUSE_OPENBLAS"] + self.extra_compile_args += ["/DUSE_OPENBLAS"] print(sys.version) loc = dirname(sys.executable) @@ -93,7 +93,7 @@ def conda_setup_for_unix(self): self.include_dirs += self.check_valid_path([join(loc, "include")]) self.library_dirs += self.check_valid_path([join(loc, "lib")]) - self.extra_compile_args = ["-O3", "-ffast-math", "-DUSE_OPENBLAS"] + self.extra_compile_args += ["-O3", "-ffast-math", "-DUSE_OPENBLAS"] self.libraries += ["fftw3", "openblas"] def on_exit_message(self, blas_lib, fftw_lib): @@ -132,8 +132,8 @@ class WindowsSetup(Setup): def __init__(self): super().__init__() - self.extra_link_args += ["-Wl"] - self.extra_compile_args = ["-DFFTW_DLL"] + self.extra_link_args += [] + self.extra_compile_args += ["-DFFTW_DLL", "/permissive-"] # if use_mkl: # self.mkl_blas_info() @@ -148,6 +148,11 @@ def __init__(self): "-O3", "-ffast-math", "-fcommon", + "-Wall", + "-Wextra", + "-Wconversion", + # "-Werror", + # "-pedantic", # "-msse4.2", # "-ftree-vectorize", # "-fopt-info-vec-all", @@ -185,6 +190,9 @@ def __init__(self): # "-Rpass-analysis=loop-vectorize", "-fvectorize", "-fcommon", + "-Wall", + "-Wextra", + "-Wconversion", ] self.extra_link_args += ["-lm"] diff --git a/src/c_lib/base/base_model.pxd b/src/c_lib/base/base_model.pxd index 797921e76..b452fd7e1 100644 --- a/src/c_lib/base/base_model.pxd +++ b/src/c_lib/base/base_model.pxd @@ -26,44 +26,44 @@ cdef extern from "angular_momentum/wigner_matrix.h": cdef extern from "schemes.h": ctypedef struct MRS_averaging_scheme: - unsigned int total_orientations + int total_orientations ctypedef struct MRS_fftw_scheme: pass MRS_averaging_scheme *MRS_create_averaging_scheme( - unsigned int integration_density, + int integration_density, bool_t allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, + int n_gamma, + int integration_volume, bool_t interpolation) MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( double *alpha, double *beta, double *weight, - unsigned int n_angles, + int n_angles, bool_t allow_4th_rank, - unsigned int n_gamma, - unsigned int position_size, + int n_gamma, + int position_size, int32_t *positions, bool_t interpolation) void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) - MRS_fftw_scheme *create_fftw_scheme(unsigned int total_orientations, - unsigned int number_of_sidebands) + MRS_fftw_scheme *create_fftw_scheme(int total_orientations, + int number_of_sidebands) void MRS_free_fftw_scheme(MRS_fftw_scheme *fftw_scheme) cdef extern from "mrsimulator.h": ctypedef struct MRS_plan: - unsigned int number_of_sidebands + int number_of_sidebands double rotor_frequency_in_Hz double rotor_angle_in_rad - MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, + MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, - double rotor_angle_in_rad, double increment, + double rotor_angle_in_rad, bool_t allow_4th_rank) void MRS_free_plan(MRS_plan *plan) void MRS_get_amplitudes_from_plan(MRS_plan *plan, bool_t refresh) @@ -110,7 +110,7 @@ cdef extern from "method.h": double increment # Increment of coordinates along the dimension. double coordinates_offset # Start coordinate of the dimension. MRS_event *events # Holds a list of events. - unsigned int n_events # The number of events. + int n_events # The number of events. MRS_dimension *MRS_create_dimensions( MRS_averaging_scheme *scheme, @@ -124,8 +124,8 @@ cdef extern from "method.h": double *rotor_frequency_in_Hz, double *rotor_angle_in_rad, int *n_events, - unsigned int n_dim, - unsigned int *number_of_sidebands) + int n_dim, + int *number_of_sidebands) void MRS_free_dimension(MRS_dimension *dimensions, int n) @@ -134,9 +134,9 @@ cdef extern from "simulation.h": void mrsimulator_core( # spectrum information and related amplitude double *spec, - double spectral_start, - double spectral_increment, - int number_of_points, + # double spectral_start, + # double spectral_increment, + # int number_of_points, site_struct *sites, coupling_struct *couplings, @@ -147,14 +147,14 @@ cdef extern from "simulation.h": int quad_second_order, # Quad theory for second order, # spin rate, spin angle and number spinning sidebands - unsigned int number_of_sidebands, + int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, float *transition_pathway, # Pointer to a list of transitions. int integration_density, - unsigned int integration_volume, # 0-octant, 1-hemisphere, 2-sphere - unsigned int interpolate_type, + int integration_volume, # 0-octant, 1-hemisphere, 2-sphere + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix, ) @@ -170,7 +170,7 @@ cdef extern from "simulation.h": MRS_dimension *dimensions, # the dimensions within method. MRS_fftw_scheme *fftw_scheme, # the fftw scheme MRS_averaging_scheme *scheme, # the powder averaging scheme - unsigned int interpolate_type, + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix, ) diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index 3762303ea..ee1b34160 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -19,12 +19,12 @@ def core_simulator(method, list transition_pathways, # Same length as spin_systems list and list transition_weights, # elements coorespond to each system int verbose=0, # for debug purpose only. - unsigned int number_of_sidebands=90, - unsigned int integration_density=72, - unsigned int decompose_spectrum=0, - unsigned int integration_volume=1, - unsigned int isotropic_interpolation=0, - unsigned int number_of_gamma_angles=1, + int number_of_sidebands=90, + int integration_density=72, + int decompose_spectrum=0, + int integration_volume=1, + int isotropic_interpolation=0, + int number_of_gamma_angles=1, bool_t interpolation=True, bool_t auto_switch=True, debug=False, @@ -56,15 +56,15 @@ def core_simulator(method, # create averaging scheme _____________________________________________________ cdef clib.MRS_averaging_scheme *averaging_scheme - cdef unsigned int position_size = 0 + cdef int position_size = 0 if user_defined: if interpolation: - position_size = np.uint32(positions.size / 3) if positions is not None else 0 + position_size = np.intc(positions.size / 3) if positions is not None else 0 else: positions = None averaging_scheme = clib.MRS_create_averaging_scheme_from_alpha_beta( - alpha=&alpha[0], beta=&beta[0], weight=&weight[0], n_angles=alpha.size, + alpha=&alpha[0], beta=&beta[0], weight=&weight[0], n_angles=np.intc(alpha.size), allow_4th_rank=allow_4th_rank, n_gamma=number_of_gamma_angles, position_size=position_size, positions=&positions[0], interpolation=interpolation ) @@ -76,7 +76,7 @@ def core_simulator(method, ) # create C spectral dimensions ________________________________________________ - cdef int n_dimension = len(method.spectral_dimensions) + cdef int n_dimension = int(len(method.spectral_dimensions)) # Define and allocate C numpy arrays n_points = 1 @@ -91,7 +91,7 @@ def core_simulator(method, cdef ndarray[int] cnt cdef ndarray[double] coord_off cdef ndarray[double] incre - cdef ndarray[unsigned int] n_dim_sidebands + cdef ndarray[int] n_dim_sidebands # Loop through dimensions and grab attributes/values in python freq_contrib = np.asarray([]) @@ -155,11 +155,11 @@ def core_simulator(method, magnetic_flux_density_in_T = np.asarray(Bo, dtype=np.float64) srfiH = np.asarray(vr, dtype=np.float64) rair = np.asarray(th, dtype=np.float64) - cnt = np.asarray(count, dtype=np.int32) + cnt = np.asarray(count, dtype=np.intc) incre = np.asarray(increment, dtype=np.float64) coord_off = np.asarray(coordinates_offset, dtype=np.float64) - n_event = np.asarray(event_i, dtype=np.int32) - n_dim_sidebands = np.asarray(dim_sidebands, dtype=np.uint32) + n_event = np.asarray(event_i, dtype=np.intc) + n_dim_sidebands = np.asarray(dim_sidebands, dtype=np.intc) # # special 1D case with 1 event. # if np.all(srfiH == 1e-3) and np.all(rair - rair[0] == 0): @@ -182,7 +182,7 @@ def core_simulator(method, norm = np.abs(np.prod(incre)) # create fftw scheme __________________________________________________________ - cdef unsigned int max_sidebands = n_dim_sidebands.max() + cdef int max_sidebands = n_dim_sidebands.max() cdef clib.MRS_fftw_scheme *fftw_scheme fftw_scheme = clib.create_fftw_scheme(averaging_scheme.total_orientations, max_sidebands) # _____________________________________________________________________________ @@ -206,7 +206,7 @@ def core_simulator(method, affine_matrix_c[3] -= affine_matrix_c[1]*affine_matrix_c[2] # sites _______________________________________________________________________________ - cdef int number_of_sites, number_of_couplings + cdef int number_of_sites, number_of_couplings, _site_idx_, _coup_idx_ cdef ndarray[int] spin_index_ij cdef ndarray[float] spin_i cdef ndarray[double] gyromagnetic_ratio_i @@ -258,7 +258,7 @@ def core_simulator(method, # sub_sites = [site for site in spin_sys.sites if site.isotope.symbol == isotope] # index_.append(index) - number_of_sites = len(spin_sys.sites) + number_of_sites = int(len(spin_sys.sites)) # ------------------------------------------------------------------------ # Site specification @@ -280,23 +280,23 @@ def core_simulator(method, # Extract and assign site information from Site objects to C structure # --------------------------------------------------------------------- - for i in range(number_of_sites): - site = spin_sys.sites[i] - spin_i[i] = site.isotope.spin - gyromagnetic_ratio_i[i] = site.isotope.gyromagnetic_ratio - one_minus_sigma_iso_ref_i[i] = site.isotope.ref_larmor_ratio - i3 = 3*i + for _site_idx_ in range(number_of_sites): + site = spin_sys.sites[_site_idx_] + spin_i[_site_idx_] = np.float32(site.isotope.spin) + gyromagnetic_ratio_i[_site_idx_] = site.isotope.gyromagnetic_ratio + one_minus_sigma_iso_ref_i[_site_idx_] = site.isotope.ref_larmor_ratio + i3 = 3 * _site_idx_ # CSA tensor if site.isotropic_chemical_shift is not None: - iso_n[i] = site.isotropic_chemical_shift + iso_n[_site_idx_] = site.isotropic_chemical_shift shielding = site.shielding_symmetric if shielding is not None: if shielding.zeta is not None: - zeta_n[i] = shielding.zeta + zeta_n[_site_idx_] = shielding.zeta if shielding.eta is not None: - eta_n[i] = shielding.eta + eta_n[_site_idx_] = shielding.eta if shielding.alpha is not None: ori_n[i3] = shielding.alpha if shielding.beta is not None: @@ -314,9 +314,9 @@ def core_simulator(method, quad = site.quadrupolar if quad is not None: if quad.Cq is not None: - Cq_e[i] = quad.Cq + Cq_e[_site_idx_] = quad.Cq if quad.eta is not None: - eta_e[i] = quad.eta + eta_e[_site_idx_] = quad.eta if quad.alpha is not None: ori_e[i3] = quad.alpha if quad.beta is not None: @@ -327,7 +327,7 @@ def core_simulator(method, if debug: print(f'Quadrupolar coupling constant (Cq) = {Cq_e[i]/1e6} MHz') print(f'Quadrupolar asymmetry (η) = {eta_e}') - print(f'Quadrupolar orientation (alpha/beta/gamma) = {ori_e}]') + print(f'Quadrupolar orientation (alpha/beta/gamma) = {ori_e}') # sites packed as c struct sites_c.number_of_sites = number_of_sites @@ -349,8 +349,8 @@ def core_simulator(method, # J-coupling couplings_c.number_of_couplings = 0 if spin_sys.couplings is not None: - number_of_couplings = len(spin_sys.couplings) - spin_index_ij = np.zeros(2*number_of_couplings, dtype=np.int32) + number_of_couplings = int(len(spin_sys.couplings)) + spin_index_ij = np.zeros(2*number_of_couplings, dtype=np.intc) iso_j = np.zeros(number_of_couplings, dtype=np.float64) zeta_j = np.zeros(number_of_couplings, dtype=np.float64) @@ -363,21 +363,21 @@ def core_simulator(method, ori_d = np.zeros(3*number_of_couplings, dtype=np.float64) # Extract and assign coupling information from Site objects to C structure - for i in range(number_of_couplings): - coupling = spin_sys.couplings[i] - spin_index_ij[2*i: 2*i+2] = coupling.site_index - i3 = 3*i + for _coup_idx_ in range(number_of_couplings): + coupling = spin_sys.couplings[_coup_idx_] + spin_index_ij[2*_coup_idx_: 2*_coup_idx_+2] = coupling.site_index + i3 = 3 * _coup_idx_ # J tensor if coupling.isotropic_j is not None: - iso_j[i] = coupling.isotropic_j + iso_j[_coup_idx_] = coupling.isotropic_j J_sym = coupling.j_symmetric if J_sym is not None: if J_sym.zeta is not None: - zeta_j[i] = J_sym.zeta + zeta_j[_coup_idx_] = J_sym.zeta if J_sym.eta is not None: - eta_j[i] = J_sym.eta + eta_j[_coup_idx_] = J_sym.eta if J_sym.alpha is not None: ori_j[i3] = J_sym.alpha if J_sym.beta is not None: @@ -389,9 +389,9 @@ def core_simulator(method, dipolar = coupling.dipolar if dipolar is not None: if dipolar.D is not None: - D_d[i] = dipolar.D + D_d[_coup_idx_] = dipolar.D if dipolar.eta is not None: - eta_d[i] = dipolar.eta + eta_d[_coup_idx_] = dipolar.eta if dipolar.alpha is not None: ori_d[i3] = dipolar.alpha if dipolar.beta is not None: @@ -399,17 +399,17 @@ def core_simulator(method, if dipolar.gamma is not None: ori_d[i3+2] = dipolar.gamma - # if debug: - # print(f'N couplings = {number_of_couplings}') - # print(f'site index J = {spin_index_ij}') - # print(f'Isotropic J = {iso_j} Hz') - # print(f'J anisotropy = {zeta_j} Hz') - # print(f'J asymmetry = {eta_j}') - # print(f'J orientation = {ori_j}') + if debug: + print(f'N couplings = {number_of_couplings}') + print(f'site index J = {spin_index_ij}') + print(f'Isotropic J = {iso_j} Hz') + print(f'J anisotropy = {zeta_j} Hz') + print(f'J asymmetry = {eta_j}') + print(f'J orientation = {ori_j}') - # print(f'Dipolar coupling constant = {D_d} Hz') - # print(f'Dipolar asymmetry = {eta_d}') - # print(f'Dipolar orientation = {ori_d}') + print(f'Dipolar coupling constant = {D_d} Hz') + print(f'Dipolar asymmetry = {eta_d}') + print(f'Dipolar orientation = {ori_d}') # couplings packed as c struct couplings_c.number_of_couplings = number_of_couplings @@ -486,6 +486,7 @@ def core_simulator(method, amp1.shape = method.shape() if gyromagnetic_ratio < 0: amp1 = np.fft.fftn(np.fft.ifftn(amp1).conj()) + amp1 = [amp1] clib.MRS_free_dimension(dimensions, n_dimension) clib.MRS_free_averaging_scheme(averaging_scheme) @@ -497,7 +498,7 @@ def core_simulator(method, @cython.boundscheck(False) @cython.wraparound(False) def get_zeeman_states(sys): - cdef int i, j, n_site = len(sys.sites) + cdef int i, j, n_site = int(len(sys.sites)) two_Ip1 = [int(2 * site.isotope.spin + 1) for site in sys.sites] spin_quantum_numbers = [ @@ -597,7 +598,7 @@ def calculate_transition_connect_weight( Return: A complex amplitude. """ cdef ndarray[double] factor = np.asarray([1, 0], dtype=np.float64) - cdef int i, n_sites = spin.size + cdef int i, n_sites = int(spin.size) cdef float m1_f, m1_i, m2_f, m2_i for i in range(n_sites): m1_f = trans1[1][i] # starting transition final state @@ -621,7 +622,7 @@ def get_Haeberlen_components( double eta, double rho): """Return random extended czjzek tensors in Haeberlen convention""" - cdef int n = expr_base_p.shape[1] + cdef int n = int(expr_base_p.shape[1]) cdef ndarray[double, ndim=2] param = np.empty((n, 2), dtype=float) clib.vm_haeberlen_components( n, &expr_base_p[0, 0], &expr_base_q[0, 0], zeta, eta, rho, ¶m[0, 0] diff --git a/src/c_lib/clib/clib.pyx b/src/c_lib/clib/clib.pyx index e42ef705f..8aeec6842 100644 --- a/src/c_lib/clib/clib.pyx +++ b/src/c_lib/clib/clib.pyx @@ -20,7 +20,7 @@ def histogram1d( ndarray[double, ndim=1] weights, ): cdef ndarray[double] hist = np.zeros(x_count, dtype=np.float64) - cdef int sample_count = sample_x.size + cdef int sample_count = int(sample_x.size) cdef int c_int64 = ctypes.sizeof(ctypes.c_int64) cdef int stride_x = sample_x.ctypes.strides[0] / c_int64 cdef int stride_w = weights.ctypes.strides[0] / c_int64 @@ -49,7 +49,7 @@ def histogram2d( int interp=0, ): cdef ndarray[double, ndim=2] hist = np.zeros((x_count, y_count), dtype=np.float64) - cdef int sample_count = sample_x.size + cdef int sample_count = int(sample_x.size) cdef int c_int64 = ctypes.sizeof(ctypes.c_int64) cdef int stride_x = sample_x.ctypes.strides[0] / c_int64 cdef int stride_y = sample_y.ctypes.strides[0] / c_int64 diff --git a/src/c_lib/include/angular_momentum/wigner_matrix.h b/src/c_lib/include/angular_momentum/wigner_matrix.h index 0372d52c2..4d4234c72 100644 --- a/src/c_lib/include/angular_momentum/wigner_matrix.h +++ b/src/c_lib/include/angular_momentum/wigner_matrix.h @@ -121,8 +121,7 @@ extern void single_wigner_rotation(const int l, const double *euler_angles, * fourth-rank wigner matrices. The length of w4 is `octant_orientations x * n_octants x 9` with 9 as the leading dimension. */ -extern void __batch_wigner_rotation(const unsigned int octant_orientations, - const unsigned int n_octants, +extern void __batch_wigner_rotation(const int octant_orientations, const int n_octants, double *wigner_2j_matrices, complex128 *R2, double *wigner_4j_matrices, complex128 *R4, complex128 *exp_Im_alpha, complex128 *w2, @@ -133,5 +132,5 @@ extern void __batch_wigner_rotation(const unsigned int octant_orientations, * The function accepts cos_angle = cos(angle). * The result is stored in exp_Im_angle as m x n matrix where m = [-4,-3,-2,-1] */ -extern void get_exp_Im_angle(const unsigned int n, const bool allow_4th_rank, - void *exp_Im_angle, double delta_alpha); +extern void get_exp_Im_angle(const int n, const bool allow_4th_rank, void *exp_Im_angle, + double delta_alpha); diff --git a/src/c_lib/include/config.h b/src/c_lib/include/config.h index cbdf45a2c..18f9ab54d 100644 --- a/src/c_lib/include/config.h +++ b/src/c_lib/include/config.h @@ -60,12 +60,13 @@ typedef float complex64[2]; // ---------------------------------------------------------------------------- // // OS base definitions -------------------------------------------------------- // -#define __int64_ long // for both mac-os and linux (unix system) - #ifdef _WIN32 // windows 32-bit or 64-bit #define __int64_ long long // #if _MSC_VER && !__INTEL_COMPILER // windows msvc compiler // #endif // end windows 32-bit or 64-bit +#else +// for both mac-os and linux (unix system) +#define __int64_ long #endif // end windows msvc compiler // ---------------------------------------------------------------------------- // diff --git a/src/c_lib/include/frequency_averaging.h b/src/c_lib/include/frequency_averaging.h index 61b7bc335..7bd85deb7 100644 --- a/src/c_lib/include/frequency_averaging.h +++ b/src/c_lib/include/frequency_averaging.h @@ -10,9 +10,8 @@ #include "simulation.h" void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, unsigned int iso_intrp, - complex128 *exp_I_phase); + double *spec, int iso_intrp, complex128 *exp_I_phase); void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, double *affine_matrix, - unsigned int iso_intrp, complex128 *exp_I_phase); + double *spec, double *affine_matrix, int iso_intrp, + complex128 *exp_I_phase); diff --git a/src/c_lib/include/histogram.h b/src/c_lib/include/histogram.h index 6eb57029f..f759a216c 100644 --- a/src/c_lib/include/histogram.h +++ b/src/c_lib/include/histogram.h @@ -9,19 +9,19 @@ extern void histogram1d_c(int x_count, const double x_min, const double x_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict weights, const int stride_w); + double *restrict sample_x, const int stride_x, + double *restrict weights, const int stride_w); extern void histogram2d_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict sample_y, const int stride_y, - const double *restrict weights, const int stride_w); + double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w); extern void histogram2d_interp_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict sample_y, const int stride_y, - const double *restrict weights, const int stride_w); + double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w); diff --git a/src/c_lib/include/interpolation.h b/src/c_lib/include/interpolation.h index 1f8d2cc94..96fd86891 100644 --- a/src/c_lib/include/interpolation.h +++ b/src/c_lib/include/interpolation.h @@ -12,96 +12,89 @@ /** * @brief Create a triangle with coordinates (f1, f2, f2) onto a 1D grid. * - * @param f1 A pointer to the coordinate f11. - * @param f2 A pointer to the coordinate f12. - * @param f3 A pointer to the coordinate f13. - * @param amp A pointer to the area of the vector. + * @param f1 Coordinate f11. + * @param f2 Coordinate f12. + * @param f3 Coordinate f13. + * @param amp Area of triangle. * @param spec A pointer to the starting index of a one-dimensional array. - * @param m0 A pointer to the number of points on the 1D grid. + * @param m0 The number of points on the 1D grid. * @param type The type of interpolation for the isotropic components. * 0. delta-interpolation, * 1. gaussian-interpolation. */ -extern void triangle_interpolation1D(double *f1, double *f2, double *f3, double *amp, - double *spec, int *m0, unsigned int iso_intrp); +extern void triangle_interpolation1D(double f1, double f2, double f3, double amp, + double *spec, int m0, int iso_intrp); -extern void triangle_interpolation1D_linear(double *f1, double *f2, double *f3, - double *amp, double *spec, int *m0); +extern void triangle_interpolation1D_linear(double f1, double f2, double f3, double amp, + double *spec, int m0); -extern void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, +extern void one_d_averaging(double *spec, const int freq_size, double *freq, double *amp_real, double *amp_imag, int dimension_count, - const unsigned int position_size, int32_t *positions, - const unsigned int nt, bool user_defined, - bool interpolation); + const int position_size, int32_t *positions, const int nt, + bool user_defined, bool interpolation); -extern void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, +extern void two_d_averaging(double *spec, const int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, - const unsigned int position_size, int32_t *positions, - int dimension0_count, int dimension1_count, - unsigned int iso_intrp, const unsigned int nt, - bool user_defined, bool interpolation); + const int position_size, int32_t *positions, + int dimension0_count, int dimension1_count, int iso_intrp, + const int nt, bool user_defined, bool interpolation); -extern void hist1d(double *spec, const unsigned int freq_size, double *freq, - double *amp, int m, const unsigned int nt); +extern void hist1d(double *spec, const int freq_size, double *freq, double *amp, int m); -extern void hist2d(double *spec, const unsigned int freq_size, double *freq_1, - double *freq_2, double *amp, int amp_stride, int m0, int m1, - const unsigned int nt); +extern void hist2d(double *spec, const int freq_size, double *freq_1, double *freq_2, + double *amp, int amp_stride, int m0, int m1); -extern void generic_2d_triangle_average(double *spec, const unsigned int freq_size, +extern void generic_2d_triangle_average(double *spec, const int freq_size, double *freq1, double *freq2, double *amp, int amp_stride, int m0, int m1, - const unsigned int position_size, - int32_t *positions, const unsigned int nt, - unsigned int iso_intrp); + const int position_size, int32_t *positions, + int iso_intrp); -extern void generic_1d_triangle_average(double *spec, const unsigned int freq_size, - double *freq, double *amp, int m, - const unsigned int position_size, - int32_t *positions, const unsigned int nt); +extern void generic_1d_triangle_average(double *spec, const int freq_size, double *freq, + double *amp, int m, const int position_size, + int32_t *positions); -extern void triangle_interpolation1D_gaussian(double *f1, double *f2, double *f3, - double *amp, double *spec, int *m0); +extern void triangle_interpolation1D_gaussian(double f1, double f2, double f3, + double amp, double *spec, int m0); /** * @brief Rasterize a vector triangle with coordinates ((f11, f21), (f12, f22), (f13, * f23)) onto a 2D grid. * - * @param f11 A pointer to the coordinate f11. - * @param f12 A pointer to the coordinate f12. - * @param f13 A pointer to the coordinate f13. - * @param f21 A pointer to the coordinate f21. - * @param f22 A pointer to the coordinate f22. - * @param f23 A pointer to the coordinate f23. - * @param amp A pointer to the area of the vector. + * @param f11 Coordinate f11. + * @param f12 Coordinate f12. + * @param f13 Coordinate f13. + * @param f21 Coordinate f21. + * @param f22 Coordinate f22. + * @param f23 Coordinate f23. + * @param amp Area of 2D triangle. * @param spec A pointer to the starting index of a two-dimensional array. * @param m0 An integer with the rows in the 2D grid. * @param m1 An integer with the columns in the 2D grid. * @param iso_intrp Linear=0 | Gaussian=1 isotropic interpolation scheme. */ -extern void triangle_interpolation2D(double *f11, double *f12, double *f13, double *f21, - double *f22, double *f23, double *amp, - double *spec, int m0, int m1, - unsigned int iso_intrp); +extern void triangle_interpolation2D(double f11, double f12, double f13, double f21, + double f22, double f23, double amp, double *spec, + int m0, int m1, int iso_intrp); /** * @brief Sum amplitudes from the triangles interpolations over the region of an octant. * The samplings over the octant is as per Alderman and Grand scheme. * * @param nt Number of triangles along the edge of the octant. - * @param freq A pointer to an array of frequencies evaluated at octant coordinates. + * @param freq Delta frequency. * @param amp A pointer to the amplitudes for the frequencies at octant coordinates. * @param stride Stride setp for the amplitudes (amp) array. * @param n_spec Number of points in the spectrum array (spec) * @param spec A pointer to the starting index of a one-dimensional array * @param iso_intrp Linear=0 | Gaussian=1 isotropic interpolation scheme. */ -void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *amp, - int stride, int n_spec, double *spec, - unsigned int iso_intrp); +void octahedronDeltaInterpolation(const int nt, double freq, double *amp, + const int stride, int n_spec, double *spec, + int iso_intrp); -extern void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, - double *amp, int stride, int m); +extern void octahedronInterpolation(double *spec, double *freq, const int nt, + double *amp, const int stride, int m); extern void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, - int nt, double *amp, int stride, int m0, int m1, - unsigned int iso_intrp); + const int nt, double *amp, const int stride, + int m0, int m1, int iso_intrp); diff --git a/src/c_lib/include/method.h b/src/c_lib/include/method.h index 17f82ef79..862d57cf8 100644 --- a/src/c_lib/include/method.h +++ b/src/c_lib/include/method.h @@ -29,7 +29,7 @@ typedef struct MRS_dimension { double increment; /**< Increment of coordinates along the dimension. */ double coordinates_offset; /**< Start coordinate of the dimension. */ MRS_event *events; /**< Holds a list of events. */ - unsigned int n_events; /**< The number of events. */ + int n_events; /**< The number of events. */ /* private attributes */ double R0_offset; // holds the isotropic offset. This is used in determining if or @@ -68,8 +68,7 @@ MRS_dimension *MRS_create_dimensions( MRS_averaging_scheme *scheme, int *count, double *coordinates_offset, double *increment, double *fractions, double *durations, unsigned char *is_spectral, double *magnetic_flux_density_in_T, double *rotor_frequency_in_Hz, - double *rotor_angle_in_rad, int *n_events, unsigned int n_dim, - unsigned int *number_of_sidebands); + double *rotor_angle_in_rad, int *n_events, int n_dim, int *number_of_sidebands); /** * @brief Free the memory allocation for the MRS dimensions. @@ -77,6 +76,6 @@ MRS_dimension *MRS_create_dimensions( * @param dimensions The pointer to an array of MRS_dimension structs. * @param n An interger defining the number of MRS_dimension structs in `dimensions`. */ -void MRS_free_dimension(MRS_dimension *dimensions, unsigned int n); +void MRS_free_dimension(MRS_dimension *dimensions, int n); #endif /* method_h */ diff --git a/src/c_lib/include/mrsimulator.h b/src/c_lib/include/mrsimulator.h index 4f536def2..286f280ac 100644 --- a/src/c_lib/include/mrsimulator.h +++ b/src/c_lib/include/mrsimulator.h @@ -34,8 +34,8 @@ */ struct MRS_plan { - unsigned int number_of_sidebands; /**< The number of sidebands to compute. */ - double rotor_frequency_in_Hz; /**< The sample rotation frequency in Hz. */ + int number_of_sidebands; /**< The number of sidebands to compute. */ + double rotor_frequency_in_Hz; /**< The sample rotation frequency in Hz. */ /** * The angle, in radians, describing the sample axis-rotation with respect to @@ -62,8 +62,8 @@ struct MRS_plan { bool copy_for_rotor_freq; // Set True if plan is copied from rotor freq update. bool allow_4th_rank; // If true, creates buffer/tables for 4th-rank tensors. bool is_static; // It true, compute static frequencies - unsigned int size; // # of angular orientations * number of sizebands. - unsigned int n_octants; // # of octants used in the orientational averaging. + int size; // number of angular orientations * number of sidebands. + int n_octants; // number of octants used in the orientational averaging. double *norm_amplitudes; // array of normalized amplitudes per orientation. double *wigner_d2m0_vector; // wigner-2j dm0 vector, n ∈ [-2, 2]. double *wigner_d4m0_vector; // wigner-4j dm0 vector, n ∈ [-4, 4]. @@ -83,15 +83,13 @@ typedef struct MRS_plan MRS_plan; * @param rotor_frequency_in_Hz The sample rotation frequency in Hz. * @param rotor_angle_in_rad The polar angle in radians with respect to the * z-axis describing the axis of rotation. - * @param increment The increment along the spectroscopic dimension in Hz. * @param allow_4th_rank When true, the plan calculates matrices for * processing the fourth-rank tensors. * @return A pointer to the MRS_plan. */ -MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, - unsigned int number_of_sidebands, +MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, - double increment, bool allow_4th_rank); + bool allow_4th_rank); /** * @brief Release the memory allocated for the given mrsimulator plan. @@ -224,7 +222,7 @@ void MRS_rotate_components_from_PAS_to_common_frame( unsigned char *freq_contrib // The pointer to freq contribs boolean. ); -extern void get_sideband_phase_components(unsigned int number_of_sidebands, +extern void get_sideband_phase_components(int number_of_sidebands, double spin_frequency, double *restrict pre_phase); diff --git a/src/c_lib/include/object_struct.h b/src/c_lib/include/object_struct.h index d8d0a4840..690bd1924 100644 --- a/src/c_lib/include/object_struct.h +++ b/src/c_lib/include/object_struct.h @@ -17,7 +17,7 @@ * Site structure is a collection of site interaction parameters. **/ struct __site_struct { - unsigned int number_of_sites; /**< Number of sites */ + int number_of_sites; /**< Number of sites */ /* Pointer to an array of spin quantum numbers for each site within a spin system. */ float *spin; @@ -71,7 +71,7 @@ typedef struct __site_struct site_struct; * Coupling structure is a collection of coupled site interaction parameters. **/ struct __coupling_struct { - unsigned int number_of_couplings; /**< Number of couplings */ + int number_of_couplings; /**< Number of couplings */ /* Pointer to an array of the site indexes for each coupled pair within a spin system, * incremented in stride of 2/pair. The array size is 2*number_of_couplings. */ diff --git a/src/c_lib/include/octahedron.h b/src/c_lib/include/octahedron.h index 683a9002e..96ea3041a 100644 --- a/src/c_lib/include/octahedron.h +++ b/src/c_lib/include/octahedron.h @@ -11,20 +11,20 @@ #include "interpolation.h" extern void octahedronGetDirectionCosineSquareAndWeightsOverOctant( - const unsigned int nt, double *restrict xr, double *restrict yr, - double *restrict zr, double *restrict amp); + const int nt, double *restrict xr, double *restrict yr, double *restrict zr, + double *restrict amp); -extern void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, +extern void octahedronGetPolarAngleTrigOverOctant(const int nt, double *restrict cos_alpha, double *restrict cos_beta, double *restrict amp); -extern void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, +extern void octahedronGetComplexExpOfPolarAngleOverOctant(const int nt, void *restrict exp_I_alpha, void *restrict exp_I_beta, double *restrict amp); -void get_total_amplitude(const unsigned int nt, double *amp, double *amp_sum); +void get_total_amplitude(const int nt, double *amp, double *amp_sum); -extern void averaging_setup(unsigned int nt, void *exp_I_alpha, void *exp_I_beta, - double *amp, bool interpolation); +extern void averaging_setup(int nt, void *exp_I_alpha, void *exp_I_beta, double *amp, + bool interpolation); diff --git a/src/c_lib/include/schemes.h b/src/c_lib/include/schemes.h index 87e72ceea..869aa5d12 100644 --- a/src/c_lib/include/schemes.h +++ b/src/c_lib/include/schemes.h @@ -36,30 +36,30 @@ * thousands of sites. */ typedef struct MRS_averaging_scheme { - unsigned int total_orientations; /**< The total number of orientations. */ + int total_orientations; /**< The total number of orientations. */ /** \privatesection */ - unsigned int integration_density; // # triangles along the edge of the octahedron. - unsigned int integration_volume; // 0-octant, 1-hemisphere, 2-sphere. - unsigned int octant_orientations; // # unique orientations on the face of an octant. - unsigned int n_gamma; // number of gamma angles - double *amplitudes; // array of amplitude scaling per orientation. - complex128 *exp_Im_alpha; // array of e^in\alpha n=[0,4] per orientation. - complex128 *exp_Im_gamma; // array of e^im\gamma m=[0,4] per orientation. - complex128 *w2; // buffer for 2nd rank frequency calculation. - complex128 *w4; // buffer for 4nd rank frequency calculation. - double *wigner_2j_matrices; // wigner-d 2j matrix per orientation. - double *wigner_4j_matrices; // wigner-d 4j matrix per orientation. - double *scratch; // scratch memory for calculations. - bool allow_4th_rank; // If true, compute wigner matrices for wigner-d 4j. - double *amps_real; // array of real amplitudes - double *amps_imag; // array of real amplitudes + int integration_density; // number of triangles along the edge of the octahedron. + int integration_volume; // 0-octant, 1-hemisphere, 2-sphere. + int octant_orientations; // number of unique orientations on the face of an octant. + int n_gamma; // number of gamma angles + double *amplitudes; // array of amplitude scaling per orientation. + complex128 *exp_Im_alpha; // array of e^in\alpha n=[0,4] per orientation. + complex128 *exp_Im_gamma; // array of e^im\gamma m=[0,4] per orientation. + complex128 *w2; // buffer for 2nd rank frequency calculation. + complex128 *w4; // buffer for 4nd rank frequency calculation. + double *wigner_2j_matrices; // wigner-d 2j matrix per orientation. + double *wigner_4j_matrices; // wigner-d 4j matrix per orientation. + double *scratch; // scratch memory for calculations. + bool allow_4th_rank; // If true, compute wigner matrices for wigner-d 4j. + double *amps_real; // array of real amplitudes + double *amps_imag; // array of real amplitudes double *phase; complex128 *exp_I_phase; - int32_t *positions; // array of triangle vertexes (faces) on mesh. - unsigned int position_size; // number of triangle vertexes (faces) on mesh - bool user_defined; // if true, the scheme is user defined - bool interpolation; // if true, use frequency triangle interpolation + int32_t *positions; // array of triangle vertexes (faces) on mesh. + int position_size; // number of triangle vertexes (faces) on mesh + bool user_defined; // if true, the scheme is user defined + bool interpolation; // if true, use frequency triangle interpolation } MRS_averaging_scheme; // typedef struct MRS_averaging_scheme; @@ -82,10 +82,9 @@ typedef struct MRS_averaging_scheme { * @param integration_volume An enumeration. 0=octant, 1=hemisphere, 2=sphere * @param interpolation If true, apply frequency triangle interpolation. */ -MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_density, - bool allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, +MRS_averaging_scheme *MRS_create_averaging_scheme(int integration_density, + bool allow_4th_rank, int n_gamma, + int integration_volume, bool interpolation); /** @@ -97,7 +96,7 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi * double. * @param weight A pointer to an array of size `n_angles` holding the weights of the * angle pair (alpha, beta) of type double. - * @param n_angles An unsigned int equal to the total number of alpha, beta, and weight + * @param n_angles An int equal to the total number of alpha, beta, and weight * values. * @param allow_4th_rank If true, the scheme also calculates matrices for fourth-rank * tensors. @@ -107,9 +106,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi * @param interpolation If true, apply frequency triangle interpolation. */ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( - double *alpha, double *beta, double *weight, unsigned int n_angles, - bool allow_4th_rank, unsigned int n_gamma, const unsigned int position_size, - int32_t *positions, bool interpolation); + double *alpha, double *beta, double *weight, int n_angles, bool allow_4th_rank, + int n_gamma, const int position_size, int32_t *positions, bool interpolation); /** * Free the memory allocated for the spatial orientation averaging scheme. @@ -132,8 +130,7 @@ typedef struct MRS_fftw_scheme { fftw_plan the_fftw_plan; // The plan for fftw routine. } MRS_fftw_scheme; -MRS_fftw_scheme *create_fftw_scheme(unsigned int total_orientations, - unsigned int number_of_sidebands); +MRS_fftw_scheme *create_fftw_scheme(int total_orientations, int number_of_sidebands); void MRS_free_fftw_scheme(MRS_fftw_scheme *fftw_scheme); diff --git a/src/c_lib/include/simulation.h b/src/c_lib/include/simulation.h index 742fcfa4c..8b8515eae 100644 --- a/src/c_lib/include/simulation.h +++ b/src/c_lib/include/simulation.h @@ -14,26 +14,26 @@ extern void mrsimulator_core( // spectrum information and related amplitude - double *spec, // Pointer to the spectrum array (complex). - double spectral_start, // The start of the frequency spectrum. - double spectral_increment, // The increment of the frequency spectrum. - int number_of_points, // Number of points on the frequency spectrum. + double *spec, // Pointer to the spectrum array (complex). + // double spectral_start, // The start of the frequency spectrum. + // double spectral_increment, // The increment of the frequency spectrum. + // int number_of_points, // Number of points on the frequency spectrum. site_struct *sites, // Pointer to sites within a spin system. coupling_struct *couplings, // Pointer to couplings within spin system. MRS_dimension *dimensions, // Pointer to the spectral dimension. int n_dimension, // Number of dimensions. int quad_second_order, // Quad theory for second order, - unsigned int number_of_sidebands, // The number of sidebands. - double rotor_frequency_in_Hz, // The rotor spin frequency. - double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis. + int number_of_sidebands, // The number of sidebands. + double rotor_frequency_in_Hz, // The rotor spin frequency. + double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis. float *transition_pathway, // Pointer to a list of transitions. double *transition_pathway_weight, // The complex weight of transition pathway. // powder orientation average int integration_density, // The number of triangle along the edge of octahedron. - unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. - unsigned int interpolate_type, unsigned char *freq_contrib, double *affine_matrix); + int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix); extern void __mrsimulator_core( // spectrum information and related amplitude @@ -51,7 +51,7 @@ extern void __mrsimulator_core( MRS_dimension *dimensions, // Pointer to MRS_dimension structure. MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. - unsigned int interpolate_type, + int interpolate_type, /** * Each event consists of the following freq contrib ordered as diff --git a/src/c_lib/include/vm.h b/src/c_lib/include/vm.h index f1719174e..6d2500501 100644 --- a/src/c_lib/include/vm.h +++ b/src/c_lib/include/vm.h @@ -20,6 +20,7 @@ static inline double absd(double a) { *((unsigned __int64_ *)&a) &= ~(1ULL << 63); return a; + // return fabs(a); } /** Arithmetic suit ======================================================== */ @@ -216,8 +217,7 @@ static inline void vm_double_square_root_inplace(int count, double *restrict x) * res = x * y */ static inline void vm_double_complex_multiply(int count, const void *restrict x, - const void *restrict y, - void *restrict res) { + const void *restrict y, void *res) { double *res_ = (double *)res; double *x_ = (double *)x; double *y_ = (double *)y; @@ -439,7 +439,7 @@ static inline void vm_double_complex_exp(int count, const void *restrict x, * res = exp(x(imag)) */ static inline void vm_double_complex_exp_imag_only(int count, const void *restrict x, - void *restrict res) { + void *res) { // x = __builtin_assume_aligned(x, 32); // res = __builtin_assume_aligned(res, 32); double *x_ = (double *)x; diff --git a/src/c_lib/include/vm_common.h b/src/c_lib/include/vm_common.h index 4a90f0383..ba7ae7c19 100644 --- a/src/c_lib/include/vm_common.h +++ b/src/c_lib/include/vm_common.h @@ -107,7 +107,7 @@ static inline void vm_double_ones(int count, double *restrict res) { * @returns values A pointer to the fft output order vector of size @p n. */ static inline double *get_FFT_order_freq(int n, double increment) { - double *vr_freq = malloc_double(n); + double *vr_freq = malloc_double((size_t)n); int m, positive_limit, negative_limit; if (n % 2 == 0) { @@ -121,7 +121,7 @@ static inline double *get_FFT_order_freq(int n, double increment) { for (m = 0; m <= positive_limit; m++) *vr_freq++ = (double)m * increment; for (m = negative_limit; m < 0; m++) *vr_freq++ = (double)m * increment; return vr_freq - n; -}; +} /** * @brief Return a vector ordered according to the fft output order. @@ -130,7 +130,7 @@ static inline double *get_FFT_order_freq(int n, double increment) { * @returns values A pointer to the fft output order vector of size @p n. */ static inline int *get_FFT_index(int n) { - int *freq_index = malloc_int(n); + int *freq_index = malloc_int((size_t)n); int m, positive_limit, negative_limit; if (n % 2 == 0) { @@ -144,4 +144,4 @@ static inline int *get_FFT_index(int n) { for (m = 0; m <= positive_limit; m++) *freq_index++ = m; for (m = negative_limit; m < 0; m++) *freq_index++ = m; return freq_index - n; -}; +} diff --git a/src/c_lib/lib/angular_momentum/wigner_element.c b/src/c_lib/lib/angular_momentum/wigner_element.c index 63a48eba5..df3e106aa 100644 --- a/src/c_lib/lib/angular_momentum/wigner_element.c +++ b/src/c_lib/lib/angular_momentum/wigner_element.c @@ -37,16 +37,17 @@ static inline double __generic_wigner_d_element(const float l, const float m1, double cx = cos(beta / 2.); double sum = 0.0, k1, k2, k3, x, y; int sign = 1; - int k; + int k, n1, n2; + float two = 2.0; for (k = 0; k <= (int)(l - m1); k++) { - k1 = (int)(l - m1 - k); - k2 = (int)(l + m2 - k); - k3 = (int)(k + m1 - m2); + k1 = (int)(l - m1 - (float)k); + k2 = (int)(l + m2 - (float)k); + k3 = (int)((float)k + m1 - m2); if (k1 >= 0 && k2 >= 0 && k3 >= 0) { - int n1 = (int)(2 * l + m2 - m1 - 2 * k); - int n2 = (int)(m1 - m2 + 2 * k); + n1 = (int)(two * l + m2 - m1 - two * (float)k); + n2 = (int)(m1 - m2 + two * (float)k); x = my_power(cx, n1); y = my_power(sx, n2); sum += sign * x * y / @@ -54,7 +55,8 @@ static inline double __generic_wigner_d_element(const float l, const float m1, } sign = -sign; } - double f = fac(l + m2) * fac(l - m2) * fac(l + m1) * fac(l - m1); + double f = fac((double)(l + m2)) * fac((double)(l - m2)) * fac((double)(l + m1)) * + fac((double)(l - m1)); sign = ((int)(m1 - m2) % 2 == 0) ? 1 : -1; f = sqrt(f); return (sign * sum * f); @@ -358,7 +360,7 @@ void transition_connect_factor(const float l, const float m1_f, const float m1_i // const double phi, double *restrict factor, int // n_sites) { // double m1_f, m1_i, m2_f, m2_i; -// complex128 *weight = malloc_complex128(n_sites); +// complex128 *weight = malloc_complex128((size_t)n_sites); // for (int i = 1; i < n_sites; i++) { // m1_i = transition_inital[i]; // m1_f = (transition_inital + n_sites)[i]; diff --git a/src/c_lib/lib/angular_momentum/wigner_matrix.c b/src/c_lib/lib/angular_momentum/wigner_matrix.c index 3f1eae6c4..3f5f044b9 100644 --- a/src/c_lib/lib/angular_momentum/wigner_matrix.c +++ b/src/c_lib/lib/angular_momentum/wigner_matrix.c @@ -17,7 +17,7 @@ complex128 NEGATIVE_IOTA = {0.0, -1.0}; // ✅ .. note: (wigner_d_matrices) tested with pytest // ......................... void wigner_d_matrices(const int l, const int n, const double *beta, double *wigner) { - complex128 *exp_I_beta = malloc_complex128(n); + complex128 *exp_I_beta = malloc_complex128((size_t)n); vm_cosine_I_sine(n, beta, exp_I_beta); wigner_d_matrices_from_exp_I_beta(l, n, false, exp_I_beta, wigner); free(exp_I_beta); @@ -741,12 +741,12 @@ void single_wigner_rotation(const int l, const double *euler_angles, const void * rotation with fourth rank wigner matrices. The length of w4 is * `octant_orientations x n_octants x 9` with 9 as the leading dimension. */ -void __batch_wigner_rotation(const unsigned int octant_orientations, - const unsigned int n_octants, double *wigner_2j_matrices, - complex128 *R2, double *wigner_4j_matrices, complex128 *R4, +void __batch_wigner_rotation(const int octant_orientations, const int n_octants, + double *wigner_2j_matrices, complex128 *R2, + double *wigner_4j_matrices, complex128 *R4, complex128 *exp_Im_alpha, complex128 *w2, complex128 *w4) { - unsigned int j, max_iter, wigner_2j_inc, wigner_4j_inc = 0; - unsigned int w2_increment, w4_increment = 0, alpha_inc; + int j, max_iter, wigner_2j_inc, wigner_4j_inc = 0; + int w2_increment, w4_increment = 0, alpha_inc; double *exp_Im_alpha_ = (double *)exp_Im_alpha; w2_increment = 3 * octant_orientations; @@ -841,13 +841,13 @@ void __batch_wigner_rotation(const unsigned int octant_orientations, * The result is stored in exp_Im_angle as m x n matrix, where m = [-4,-3,-2,-1] * Since only negative m's are computed, the following computes exp(Im angle) */ -void get_exp_Im_angle(const unsigned int n, const bool allow_4th_rank, - void *exp_Im_angle, double delta_alpha) { +void get_exp_Im_angle(const int n, const bool allow_4th_rank, void *exp_Im_angle, + double delta_alpha) { double *exp_Im_angle_ = (double *)exp_Im_angle; // The complex array is interpreted as alternating real and imag double array. // The index s_i = i*n of complex array is now at index 2*i*n. - unsigned int m_3 = 2 * n, m_2 = 4 * n, m_1 = 6 * n; // m_4 is 0 + int m_3 = 2 * n, m_2 = 4 * n, m_1 = 6 * n; // m_4 is 0 // exp(-2 I angle) vm_double_complex_multiply(n, &exp_Im_angle_[m_1], &exp_Im_angle_[m_1], @@ -863,8 +863,8 @@ void get_exp_Im_angle(const unsigned int n, const bool allow_4th_rank, } if (delta_alpha != 0.0) { - double *exp_da = malloc_double(8); - unsigned int m_4p = 8 * n, m_3p = 10 * n, m_2p = 12 * n, m_1p = 14 * n; + double *exp_da = malloc_double((size_t)8); + int m_4p = 8 * n, m_3p = 10 * n, m_2p = 12 * n, m_1p = 14 * n; exp_da[0] = cos(delta_alpha); exp_da[1] = sin(delta_alpha); diff --git a/src/c_lib/lib/frequency_averaging.c b/src/c_lib/lib/frequency_averaging.c index fc2c5557f..8c9b41906 100644 --- a/src/c_lib/lib/frequency_averaging.c +++ b/src/c_lib/lib/frequency_averaging.c @@ -64,10 +64,9 @@ static inline void sideband_amplitude(int npts, int n_octant, complex128 *a11, } void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, unsigned int iso_intrp, - complex128 *exp_I_phase) { - unsigned int i, j, k1, address, ptr, gamma_idx; - unsigned int nt = scheme->integration_density, npts = scheme->octant_orientations; + double *spec, int iso_intrp, complex128 *exp_I_phase) { + int i, j, k1, address, ptr, gamma_idx; + int nt = scheme->integration_density, npts = scheme->octant_orientations; double offset_0, offset; double *freq, *phase_ptr, *amps = dimensions->freq_amplitude; @@ -120,9 +119,9 @@ void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * vm_double_multiply(scheme->total_orientations, phase_ptr + 1, 2, &s[k1], amps_imag); while (j++ < planA->n_octants) { - octahedronDeltaInterpolation(nt, &offset, &s_real[address], 1, + octahedronDeltaInterpolation(nt, offset, &s_real[address], 1, dimensions->count, spec, iso_intrp); - octahedronDeltaInterpolation(nt, &offset, &s_imag[address], 1, + octahedronDeltaInterpolation(nt, offset, &s_imag[address], 1, dimensions->count, spec + 1, iso_intrp); address += npts; } @@ -158,10 +157,10 @@ void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * } void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - double *spec, double *affine_matrix, - unsigned int iso_intrp, complex128 *exp_I_phase) { - unsigned int i, k, j, address, gamma_idx; - unsigned int npts = scheme->octant_orientations, ptr; + double *spec, double *affine_matrix, int iso_intrp, + complex128 *exp_I_phase) { + int i, k, j, address, gamma_idx; + int npts = scheme->octant_orientations, ptr; int *fft1_index, *fft2_index, n_min, n_max; MRS_plan *planA, *planB; @@ -170,7 +169,7 @@ void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * double offset0, offset1, offsetA, offsetB; double *freq0, *freq1; double norm0, norm1; - complex128 *res_t = malloc_complex128(npts); + complex128 *res_t = malloc_complex128((size_t)npts); bool user_defined = scheme->user_defined, interpolation = scheme->interpolation; offset0 = dimensions[0].R0_offset; diff --git a/src/c_lib/lib/histogram.c b/src/c_lib/lib/histogram.c index 380c8be54..9cefc1e9e 100644 --- a/src/c_lib/lib/histogram.c +++ b/src/c_lib/lib/histogram.c @@ -13,9 +13,8 @@ #include void histogram1d_c(int x_count, const double x_min, const double x_max, - double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict weights, const int stride_w) { + double *restrict hist, int sample_count, double *restrict sample_x, + const int stride_x, double *restrict weights, const int stride_w) { int ix, counter = sample_count; double temp_sample, normx; double *x_pt = sample_x, *w_pt = weights; @@ -35,10 +34,9 @@ void histogram1d_c(int x_count, const double x_min, const double x_max, void histogram2d_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, - int sample_count, const double *restrict sample_x, - const int stride_x, const double *restrict sample_y, - const int stride_y, const double *restrict weights, - const int stride_w) { + int sample_count, double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w) { int ix, iy, hist_index, counter = sample_count; double temp_sample_x, temp_sample_y, normx, normy; double *x_pt = sample_x, *y_pt = sample_y, *w_pt = weights; @@ -66,9 +64,9 @@ void histogram2d_c(int x_count, const double x_min, const double x_max, int y_co void histogram2d_interp_c(int x_count, const double x_min, const double x_max, int y_count, const double y_min, const double y_max, double *restrict hist, int sample_count, - const double *restrict sample_x, const int stride_x, - const double *restrict sample_y, const int stride_y, - const double *restrict weights, const int stride_w) { + double *restrict sample_x, const int stride_x, + double *restrict sample_y, const int stride_y, + double *restrict weights, const int stride_w) { int ix, iy, p, counter = sample_count; bool left; double temp_sample_x, temp_sample_y, normx, normy; diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index 633fd078a..78f9f14de 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -23,59 +23,59 @@ * @param amp The area corresponding to the frequency coordinate. * @param spec1 A pointer to the intensity vector. */ -static void inline delta_fn_linear_interpolation(const double *freq, const int *points, - double *amp, double *spec1) { - int p = (int)*freq; - if (p >= *points || p < 0) return; +inline static void delta_fn_linear_interpolation(const double freq, const int points, + double amp, double *spec1) { + int p = (int)freq; + if (p >= points || p < 0) return; double diff, delta; bool left; - diff = *freq - (double)p; + diff = freq - (double)p; delta = diff - 0.5; double *spec = &spec1[2 * p]; // Do not interpolate the intensity if the difference < 1.0e-6. // Ensures that the sideband dimension frequencies do not get interpolated. if (absd(delta) < TOL) { - *spec += *amp; + *spec += amp; return; } // Linear interpolation left = delta < 0; if (left) { - if (p != 0) *(spec - 2) -= *amp * delta; - *spec += *amp * (1.0 + delta); + if (p != 0) *(spec - 2) -= amp * delta; + *spec += amp * (1.0 + delta); } else { - if (p + 1 != *points) *(spec + 2) += *amp * delta; - *spec += *amp * (1.0 - delta); + if (p + 1 != points) *(spec + 2) += amp * delta; + *spec += amp * (1.0 - delta); } } -static void inline delta_fn_gauss_interpolation(const double *freq, const int *points, - double *amp, double *spec) { +inline static void delta_fn_gauss_interpolation(const double freq, const int points, + double amp, double *spec) { double res, w, a0, a1, a2, a3, a4, sum, temp; - int p = (int)(floor(*freq)), index, ip2, ip1, im1, im2, pad = 2, p2 = 2 * p; - if (p >= *points + pad || p < -pad + 1) return; + int p = (int)(floor(freq)), index, ip2, ip1, im1, im2, pad = 2, p2 = 2 * p; + if (p >= points + pad || p < -pad + 1) return; // for sideband delta freq. It avoids round-off errors. - res = *freq - (double)p; - if (absd(res - 0.5) < TOL && p >= 0 && p < *points) { - spec[p2] += *amp; + res = freq - (double)p; + if (absd(res - 0.5) < TOL && p >= 0 && p < points) { + spec[p2] += amp; return; } - p = (int)(floor(*freq - 0.5)); + p = (int)(floor(freq - 0.5)); p2 = 2 * p; - res = *freq - (double)p - 0.5; + res = freq - (double)p - 0.5; res *= gauss_table_precision_inverse; index = (int)res; w = res - (double)index; - ip2 = 2 * gauss_table_precision_inverse - index; - ip1 = gauss_table_precision_inverse - index; - im1 = gauss_table_precision_inverse + index; - im2 = 2 * gauss_table_precision_inverse + index; + ip2 = (int)(2 * gauss_table_precision_inverse - index); + ip1 = (int)(gauss_table_precision_inverse - index); + im1 = (int)(gauss_table_precision_inverse + index); + im2 = (int)(2 * gauss_table_precision_inverse + index); a0 = lerp_plus(w, im2); a1 = lerp_plus(w, im1); @@ -89,18 +89,18 @@ static void inline delta_fn_gauss_interpolation(const double *freq, const int *p sum += a3; sum += a4; - temp = *amp / sum; + temp = amp / sum; // double *spec = &spec1[2 * p]; - if (p - 2 >= 0 && p - 2 < *points) spec[p2 - 4] += temp * a0; - if (p - 1 >= 0 && p - 1 < *points) spec[p2 - 2] += temp * a1; - if (p >= 0 && p < *points) spec[p2] += temp * a2; - if (p + 1 >= 0 && p + 1 < *points) spec[p2 + 2] += temp * a3; - if (p + 2 >= 0 && p + 2 < *points) spec[p2 + 4] += temp * a4; + if (p - 2 >= 0 && p - 2 < points) spec[p2 - 4] += temp * a0; + if (p - 1 >= 0 && p - 1 < points) spec[p2 - 2] += temp * a1; + if (p >= 0 && p < points) spec[p2] += temp * a2; + if (p + 1 >= 0 && p + 1 < points) spec[p2 + 2] += temp * a3; + if (p + 2 >= 0 && p + 2 < points) spec[p2 + 4] += temp * a4; } /** - * @brief Get the clipping conditions object + * @brief Get the clipping conditions object. Update the value of p, pmid, and pmax * * @param p Start index of the triangle. * @param pmid Apex index of the triangle. @@ -108,19 +108,19 @@ static void inline delta_fn_gauss_interpolation(const double *freq, const int *p * @param points Max attainable index. * @param clips Pointer to bool array to store the clip conditions. */ -static inline void get_clipping_conditions(int *p, int *pmid, int *pmax, int *points, +static inline void get_clipping_conditions(int *p, int *pmid, int *pmax, int points, bool *clips) { *clips = *p <= 0; // left triangle left clip *p = *clips++ ? 0 : *p; - *clips = *pmid > *points - 1; // left triangle right clip - *pmid = *clips++ ? *points - 1 : *pmid; + *clips = *pmid > points - 1; // left triangle right clip + *pmid = *clips++ ? points - 1 : *pmid; *clips = *pmid <= 0; // right triangle left clip *pmid = *clips++ ? 0 : *pmid; - *clips = *pmax > *points - 1; // right triangle right clip - *pmax = *clips ? *points - 1 : *pmax; + *clips = *pmax > points - 1; // right triangle right clip + *pmax = *clips ? points - 1 : *pmax; } /** @@ -210,20 +210,20 @@ static inline void right_triangle_interpolate(int p, int pmax, bool l_clip, bool *(spec += 2) += (!r_clip) ? df * df * 0.5 * df1 : diff - df1; } -static inline void __triangle_interpolation(double *freq1, double *freq2, double *freq3, - double *amp, double *spec, int *points) { +static inline void __triangle_interpolation(double freq1, double freq2, double freq3, + double amp, double *spec, int points) { int p, pmid, pmax, i, j; double top, t; - p = (int)(*freq1); + p = (int)(freq1); // check if the three points lie within a bin interval. - if (p == (int)(*freq2) && p == (int)(*freq3)) { - if (p >= *points || p < 0) return; - spec[2 * p] += *amp; + if (p == (int)(freq2) && p == (int)(freq3)) { + if (p >= points || p < 0) return; + spec[2 * p] += amp; return; } - double f[3] = {freq1[0], freq2[0], freq3[0]}; + double f[3] = {freq1, freq2, freq3}; // arrange the numbers in ascending order (sort) @@ -239,14 +239,14 @@ static inline void __triangle_interpolation(double *freq1, double *freq2, double // if min frequency is higher than the last bin, return p = (int)f[0]; - if (p >= *points) return; + if (p >= points) return; // if max frequency is lower than the first bin, return pmax = (int)f[2]; if (pmax < 0) return; pmid = (int)f[1]; - top = *amp * 2.0 / (f[2] - f[0]); + top = amp * 2.0 / (f[2] - f[0]); bool clips[4]; get_clipping_conditions(&p, &pmid, &pmax, points, clips); @@ -260,29 +260,32 @@ static inline void __triangle_interpolation(double *freq1, double *freq2, double } } -void triangle_interpolation1D(double *freq1, double *freq2, double *freq3, double *amp, - double *spec, int *points, unsigned int iso_intrp) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) { - if (iso_intrp == 0) return delta_fn_linear_interpolation(freq1, points, amp, spec); - if (iso_intrp == 1) return delta_fn_gauss_interpolation(freq1, points, amp, spec); +void triangle_interpolation1D(double freq1, double freq2, double freq3, double amp, + double *spec, int points, int iso_intrp) { + if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { + if (iso_intrp == 0) delta_fn_linear_interpolation(freq1, points, amp, spec); + if (iso_intrp == 1) delta_fn_gauss_interpolation(freq1, points, amp, spec); + } else { + __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); } - __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); } -void triangle_interpolation1D_linear(double *freq1, double *freq2, double *freq3, - double *amp, double *spec, int *points) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) - return delta_fn_linear_interpolation(freq1, points, amp, spec); - - __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); +void triangle_interpolation1D_linear(double freq1, double freq2, double freq3, + double amp, double *spec, int points) { + if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { + delta_fn_linear_interpolation(freq1, points, amp, spec); + } else { + __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); + } } -void triangle_interpolation1D_gaussian(double *freq1, double *freq2, double *freq3, - double *amp, double *spec, int *points) { - if (absd(*freq1 - *freq2) < TOL && absd(*freq1 - *freq3) < TOL) - return delta_fn_gauss_interpolation(freq1, points, amp, spec); - - __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); +void triangle_interpolation1D_gaussian(double freq1, double freq2, double freq3, + double amp, double *spec, int points) { + if (absd(freq1 - freq2) < TOL && absd(freq1 - freq3) < TOL) { + delta_fn_gauss_interpolation(freq1, points, amp, spec); + } else { + __triangle_interpolation(freq1, freq2, freq3, amp, spec, points); + } } /** @@ -352,20 +355,19 @@ void triangle_interpolation1D_gaussian(double *freq1, double *freq2, double *fre // / \ _____ \. // / \_____ \. // f10 ------------------ f11 bottom line -static inline void quadrilateral_bin(double *f00, double *f11, double *f10, double *f01, +static inline void quadrilateral_bin(double f00, double f11, double f10, double f01, double top, double bottom, double total, - double amp, double *spec, int m1, - unsigned int iso_intrp) { + double amp, double *spec, int m1, int iso_intrp) { double amp_bottom, amp_top, norm; if (total != 0) { norm = amp / total; amp_bottom = bottom * norm; amp_top = top * norm; - triangle_interpolation1D(f00, f11, f10, &_bottom, spec, &m1, iso_intrp); - triangle_interpolation1D(f00, f11, f01, &_top, spec, &m1, iso_intrp); + triangle_interpolation1D(f00, f11, f10, amp_bottom, spec, m1, iso_intrp); + triangle_interpolation1D(f00, f11, f01, amp_top, spec, m1, iso_intrp); } else { - triangle_interpolation1D(f00, f11, f10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f00, f11, f10, amp, spec, m1, iso_intrp); } } @@ -373,8 +375,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, bool r_clip, bool r_clip_r, double top, double *f1, double *f2, double *x10, double *x11, int m1, - double *spec, - unsigned int iso_intrp) { + double *spec, int iso_intrp) { double f10, slope_diff, abs_sdiff, abs_sdiff_2, temp, df1, diff; double amp, denom, line_up, line_down; double x00, x01, f01_slope, f02_slope; @@ -394,11 +395,11 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, *x11 = f02_slope * f10 + f2[0]; if (!r_clip && !l_clip) { amp = f10 * top * 0.5; - triangle_interpolation1D(&f2[0], x11, x10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f2[0], *x11, *x10, amp, spec, m1, iso_intrp); } if (!r_clip_r && l_clip) { amp = f1[1] * (f10 - f1[0]) * 0.5 * df1; - triangle_interpolation1D(&f2[0], x11, x10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f2[0], *x11, *x10, amp, spec, m1, iso_intrp); } return; } @@ -413,7 +414,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, *x11 = f02_slope * diff + f2[0]; if (!l_clip) { amp = 0.5 * diff * diff * df1; - triangle_interpolation1D(&f2[0], x11, x10, &, spec, &m1, iso_intrp); + triangle_interpolation1D(f2[0], *x11, *x10, amp, spec, m1, iso_intrp); } else { amp = (diff - 0.5) * df1; x00 = *x10 - f01_slope; @@ -421,7 +422,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, line_down = absd(*x11 - *x10); line_up = absd(x01 - x00); denom = line_down + line_up; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } p++; @@ -441,7 +442,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, x01 = *x11; *x10 += f01_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); line_up += abs_sdiff; line_down += abs_sdiff; @@ -459,7 +460,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, *x11 = f02_slope * f10 + f2[0]; line_down = absd(*x11 - *x10); denom = line_up + line_down; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } else { diff += df1; @@ -467,7 +468,7 @@ static inline void lower_triangle_interpolation_2d(int p, int pmid, bool l_clip, x01 = *x11; *x10 += f01_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); } } @@ -476,8 +477,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, bool r_clip, bool r_clip_l, double top, double *f1, double *f2, double *x10, double *x11, int m1, - double *spec, - unsigned int iso_intrp) { + double *spec, int iso_intrp) { double f21, slope_diff, abs_sdiff, abs_sdiff_2, temp, df2, diff; double amp, denom, line_up, line_down; double x00, x01, f12_slope, f02_slope; @@ -496,11 +496,11 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, if (p == pmax) { if (!r_clip) { amp = f21 * top * 0.5; - triangle_interpolation1D(x11, &f2[1], &f2[2], &, spec, &m1, iso_intrp); + triangle_interpolation1D(*x11, f2[1], f2[2], amp, spec, m1, iso_intrp); } if (r_clip && !r_clip_l) { amp = (f21 - diff) * (diff + f21) * 0.5 * df2; - triangle_interpolation1D(x11, &f2[1], &f2[2], &, spec, &m1, iso_intrp); + triangle_interpolation1D(*x11, f2[1], f2[2], amp, spec, m1, iso_intrp); } return; } @@ -518,7 +518,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, line_up = absd(x01 - x00); line_down = absd(*x11 - *x10); denom = line_down + line_up; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } else { amp = (diff + 0.5) * df2; @@ -529,7 +529,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, line_up = absd(x01 - x00); line_down = absd(*x11 - *x10); denom = line_down + line_up; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, amp, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, amp, spec, m1, iso_intrp); } p++; @@ -547,7 +547,7 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, x01 = *x11; *x10 += f12_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); line_up -= abs_sdiff; line_down -= abs_sdiff; @@ -562,65 +562,64 @@ static inline void upper_triangle_interpolation_2d(int p, int pmax, bool l_clip, amp = temp * temp * 0.5 * df2; x01 = *x11; *x11 = f2[2]; - triangle_interpolation1D(x10, x11, &x01, &, spec, &m1, iso_intrp); + triangle_interpolation1D(*x10, *x11, x01, amp, spec, m1, iso_intrp); } else { diff -= df2; x00 = *x10; x01 = *x11; *x10 += f12_slope; *x11 += f02_slope; - quadrilateral_bin(&x00, x11, x10, &x01, line_up, line_down, denom, diff, spec, m1, + quadrilateral_bin(x00, *x11, *x10, x01, line_up, line_down, denom, diff, spec, m1, iso_intrp); } } -void triangle_interpolation2D(double *freq11, double *freq12, double *freq13, - double *freq21, double *freq22, double *freq23, - double *amp, double *spec, int m0, int m1, - unsigned int iso_intrp) { +void triangle_interpolation2D(double freq11, double freq12, double freq13, + double freq21, double freq22, double freq23, double amp, + double *spec, int m0, int m1, int iso_intrp) { double top, t1, t2, diff, temp, n_i; int p, pmid, pmax, i, j; - double freq10_01, freq11_02; + double freq10_01 = 0.0, freq11_02 = 0.0; - p = (int)(*freq11); + p = (int)freq11; - if (absd(freq11[0] - freq12[0]) < TOL && absd(freq11[0] - freq13[0]) < TOL) { + if (absd(freq11 - freq12) < TOL && absd(freq11 - freq13) < TOL) { if (p >= m0 || p < 0) return; - diff = freq11[0] - (double)p; + diff = freq11 - (double)p; n_i = 0.5; if (absd(diff - n_i) < TOL) { - triangle_interpolation1D(freq21, freq22, freq23, amp, &spec[2 * p * m1], &m1, + triangle_interpolation1D(freq21, freq22, freq23, amp, &spec[2 * p * m1], m1, iso_intrp); return; } if (diff < n_i) { if (p != 0) { - temp = *amp * (n_i - diff); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * (p - 1) * m1], - &m1, iso_intrp); + temp = amp * (n_i - diff); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * (p - 1) * m1], + m1, iso_intrp); } - temp = *amp * (n_i + diff); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * p * m1], &m1, + temp = amp * (n_i + diff); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * p * m1], m1, iso_intrp); return; } if (diff > n_i) { if (p + 1 != m0) { - temp = *amp * (diff - n_i); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * (p + 1) * m1], - &m1, iso_intrp); + temp = amp * (diff - n_i); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * (p + 1) * m1], + m1, iso_intrp); } - temp = *amp * (1 + n_i - diff); - triangle_interpolation1D(freq21, freq22, freq23, &temp, &spec[2 * p * m1], &m1, + temp = amp * (1 + n_i - diff); + triangle_interpolation1D(freq21, freq22, freq23, temp, &spec[2 * p * m1], m1, iso_intrp); return; } return; } - double f1[3] = {freq11[0], freq12[0], freq13[0]}; - double f2[3] = {freq21[0], freq22[0], freq23[0]}; + double f1[3] = {freq11, freq12, freq13}; + double f2[3] = {freq21, freq22, freq23}; // arrange the numbers in ascending order for (j = 1; j <= 2; j++) { @@ -645,10 +644,10 @@ void triangle_interpolation2D(double *freq11, double *freq12, double *freq13, if (pmax < 0) return; pmid = (int)f1[1]; - top = *amp * 2.0 / (f1[2] - f1[0]); + top = amp * 2.0 / (f1[2] - f1[0]); bool clips[4]; - get_clipping_conditions(&p, &pmid, &pmax, &m0, clips); + get_clipping_conditions(&p, &pmid, &pmax, m0, clips); if (f1[1] >= 0.0) { lower_triangle_interpolation_2d(p, pmid, clips[0], clips[1], clips[2], top, f1, f2, &freq10_01, &freq11_02, m1, spec, iso_intrp); @@ -733,7 +732,7 @@ void rasterization(double *grid, double *v0, double *v1, double *v2, int rows, // for sec in self.clip(0, 1, 1, 0): // v1 = Point(sec[0], sec[1]) // v0 = Point(sec[2], sec[3]) -// if abs(v0.x - 1) < eps and abs(v1.x - 1) < eps \ +// if abs(v0.x - 1) < eps and abs(v1.x - 1) < eps // or abs(v0.y - 1) < eps and abs(v1.y - 1) < eps: // continue @@ -811,11 +810,11 @@ void rasterization(double *grid, double *v0, double *v1, double *v2, int rows, // self.x0 + t1*xdelta, self.y0 + t1*ydelta) // return [clipped_line] -void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *amp, - int stride, int n_spec, double *spec, - unsigned int iso_intrp) { +void octahedronDeltaInterpolation(const int nt, double freq, double *amp, + const int stride, int n_spec, double *spec, + int iso_intrp) { int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; - unsigned int int_i_stride = 0, int_j_stride = 0; + int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address; local_index = nt - 1; @@ -838,14 +837,14 @@ void octahedronDeltaInterpolation(const unsigned int nt, double *freq, double *a int_i_stride += stride; int_j_stride += stride; } - if (iso_intrp == 0) return delta_fn_linear_interpolation(freq, &n_spec, &1, spec); - if (iso_intrp == 1) return delta_fn_gauss_interpolation(freq, &n_spec, &1, spec); + if (iso_intrp == 0) delta_fn_linear_interpolation(freq, n_spec, amp1, spec); + if (iso_intrp == 1) delta_fn_gauss_interpolation(freq, n_spec, amp1, spec); } -void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, - double *amp, int stride, int m) { +void octahedronInterpolation(double *spec, double *freq, const int nt, double *amp, + const int stride, int m) { int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; - unsigned int int_i_stride = 0, int_j_stride = 0; + int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq_address; /* Interpolate between 1d points by setting up triangles of unit area */ @@ -859,12 +858,12 @@ void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, temp = amp[int_i_stride + stride] + amp_address[int_j_stride]; amp1 = temp + amp[int_i_stride]; - __triangle_interpolation(&freq[i], &freq[i + 1], &freq_address[j], &1, spec, &m); + __triangle_interpolation(freq[i], freq[i + 1], freq_address[j], amp1, spec, m); if (i < local_index) { temp += amp_address[int_j_stride + stride]; - __triangle_interpolation(&freq[i + 1], &freq_address[j], &freq_address[j + 1], - &temp, spec, &m); + __triangle_interpolation(freq[i + 1], freq_address[j], freq_address[j + 1], temp, + spec, m); } else { local_index = j + nt; i++; @@ -877,11 +876,9 @@ void octahedronInterpolation(double *spec, double *freq, const unsigned int nt, } } -void generic_1d_triangle_interpolation(double *spec, const unsigned int freq_size, - double *freq, double *amp, int m, - const unsigned int position_size, - int32_t *positions) { - unsigned int pos_size = position_size; +void generic_1d_triangle_interpolation(double *spec, double *freq, double *amp, int m, + const int position_size, int32_t *positions) { + int pos_size = position_size; int32_t p1, p2, p3; double amp_sum; @@ -891,61 +888,58 @@ void generic_1d_triangle_interpolation(double *spec, const unsigned int freq_siz p3 = *positions++; // we do amp_sum because amps are already scaled to account for the factor 3 amp_sum = amp[p1] + amp[p2] + amp[p3]; - __triangle_interpolation(&freq[p1], &freq[p2], &freq[p3], &_sum, spec, &m); + __triangle_interpolation(freq[p1], freq[p2], freq[p3], amp_sum, spec, m); } } -void hist1d(double *spec, const unsigned int freq_size, double *freq, double *amp, - int m, const unsigned int nt) { - unsigned int i = 0, ix; +void hist1d(double *spec, const int freq_size, double *freq, double *amp, int m) { + int i = 0, ix; double temp_freq; for (i = 0; i < freq_size; i++) { temp_freq = freq[i]; if (temp_freq >= 0 && temp_freq < m) { - ix = (unsigned int)temp_freq; + ix = (int)temp_freq; spec[2 * ix] += amp[i]; } } } -void one_d_averaging(double *spec, const unsigned int freq_size, double *freq, - double *amp_real, double *amp_imag, int dimension_count, - const unsigned int position_size, int32_t *positions, - const unsigned int nt, bool user_defined, bool interpolation) { +void one_d_averaging(double *spec, const int freq_size, double *freq, double *amp_real, + double *amp_imag, int dimension_count, const int position_size, + int32_t *positions, const int nt, bool user_defined, + bool interpolation) { if (!user_defined) { if (interpolation) { octahedronInterpolation(spec, freq, nt, amp_real, 1, dimension_count); octahedronInterpolation(spec + 1, freq, nt, amp_imag, 1, dimension_count); } else { - hist1d(spec, freq_size, freq, amp_real, dimension_count, nt); - hist1d(spec + 1, freq_size, freq, amp_imag, dimension_count, nt); + hist1d(spec, freq_size, freq, amp_real, dimension_count); + hist1d(spec + 1, freq_size, freq, amp_imag, dimension_count); } } else { generic_1d_triangle_average(spec, freq_size, freq, amp_real, dimension_count, - position_size, positions, nt); + position_size, positions); generic_1d_triangle_average(spec + 1, freq_size, freq, amp_imag, dimension_count, - position_size, positions, nt); + position_size, positions); } } -void generic_1d_triangle_average(double *spec, const unsigned int freq_size, - double *freq, double *amp, int m, - const unsigned int position_size, int32_t *positions, - const unsigned int nt) { +void generic_1d_triangle_average(double *spec, const int freq_size, double *freq, + double *amp, int m, const int position_size, + int32_t *positions) { if (positions == NULL) { - hist1d(spec, freq_size, freq, amp, m, nt); + hist1d(spec, freq_size, freq, amp, m); } else { - generic_1d_triangle_interpolation(spec, freq_size, freq, amp, m, position_size, - positions); + generic_1d_triangle_interpolation(spec, freq, amp, m, position_size, positions); } } -void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, - double *freq2, double *amp, int amp_stride, - const unsigned int position_size, int32_t *positions, - int dimension0_count, int dimension1_count, unsigned int iso_intrp, - const unsigned int nt, bool user_defined, bool interpolation) { +void two_d_averaging(double *spec, const int freq_size, double *freq1, double *freq2, + double *amp, int amp_stride, const int position_size, + int32_t *positions, int dimension0_count, int dimension1_count, + int iso_intrp, const int nt, bool user_defined, + bool interpolation) { if (!user_defined) { if (interpolation) { // Perform tenting on every sideband order over all orientations @@ -953,21 +947,20 @@ void two_d_averaging(double *spec, const unsigned int freq_size, double *freq1, dimension0_count, dimension1_count, iso_intrp); } else { hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, dimension0_count, - dimension1_count, nt); + dimension1_count); } } else { generic_2d_triangle_average(spec, freq_size, freq1, freq2, amp, amp_stride, dimension0_count, dimension1_count, position_size, - positions, nt, iso_intrp); + positions, iso_intrp); } } -void generic_2d_triangle_interpolation(double *spec, const unsigned int freq_size, - double *freq1, double *freq2, double *amp, - int amp_stride, const unsigned int position_size, - int32_t *positions, int m0, int m1, - unsigned int iso_intrp) { - unsigned int pos_size = position_size; +void generic_2d_triangle_interpolation(double *spec, double *freq1, double *freq2, + double *amp, int amp_stride, + const int position_size, int32_t *positions, + int m0, int m1, int iso_intrp) { + int pos_size = position_size; int32_t p1, p2, p3; double amp_sum; @@ -977,46 +970,45 @@ void generic_2d_triangle_interpolation(double *spec, const unsigned int freq_siz p3 = *positions++; // we do amp_sum because amps are already scaled to account for the factor 3 amp_sum = amp[amp_stride * p1] + amp[amp_stride * p2] + amp[amp_stride * p3]; - triangle_interpolation2D(&freq1[p1], &freq1[p2], &freq1[p3], &freq2[p1], &freq2[p2], - &freq2[p3], &_sum, spec, m0, m1, iso_intrp); + triangle_interpolation2D(freq1[p1], freq1[p2], freq1[p3], freq2[p1], freq2[p2], + freq2[p3], amp_sum, spec, m0, m1, iso_intrp); } } -void hist2d(double *spec, const unsigned int freq_size, double *freq1, double *freq2, - double *amp, int amp_stride, int m0, int m1, const unsigned int nt) { - unsigned int i = 0, ix, iy, hist_index; +void hist2d(double *spec, const int freq_size, double *freq1, double *freq2, + double *amp, int amp_stride, int m0, int m1) { + int i = 0, ix, iy, hist_index; double temp_freq_1, temp_freq_2; for (i = 0; i < freq_size; i++) { temp_freq_1 = freq1[i]; temp_freq_2 = freq2[i]; if (temp_freq_1 >= 0 && temp_freq_1 < m0 && temp_freq_2 >= 0 && temp_freq_2 < m1) { - ix = (unsigned int)temp_freq_1; - iy = (unsigned int)temp_freq_2; + ix = (int)temp_freq_1; + iy = (int)temp_freq_2; hist_index = iy + m1 * ix; spec[2 * hist_index] += amp[i * amp_stride]; } } } -void generic_2d_triangle_average(double *spec, const unsigned int freq_size, - double *freq1, double *freq2, double *amp, - int amp_stride, int m0, int m1, - const unsigned int position_size, int32_t *positions, - const unsigned int nt, unsigned int iso_intrp) { +void generic_2d_triangle_average(double *spec, const int freq_size, double *freq1, + double *freq2, double *amp, int amp_stride, int m0, + int m1, const int position_size, int32_t *positions, + int iso_intrp) { if (positions == NULL) { - hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, m0, m1, nt); + hist2d(spec, freq_size, freq1, freq2, amp, amp_stride, m0, m1); } else { - generic_2d_triangle_interpolation(spec, freq_size, freq1, freq2, amp, amp_stride, + generic_2d_triangle_interpolation(spec, freq1, freq2, amp, amp_stride, position_size, positions, m0, m1, iso_intrp); } } -void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, int nt, - double *amp, int stride, int m0, int m1, - unsigned int iso_intrp) { +void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, const int nt, + double *amp, const int stride, int m0, int m1, + int iso_intrp) { int i = 0, j = 0, local_index, n_pts = (nt + 1) * (nt + 2) / 2; - unsigned int int_i_stride = 0, int_j_stride = 0; + int int_i_stride = 0, int_j_stride = 0; double amp1, temp, *amp_address, *freq1_address, *freq2_address; /* Interpolate between 1d points by setting up triangles of unit area */ @@ -1032,15 +1024,15 @@ void octahedronInterpolation2D(double *spec, double *freq1, double *freq2, int n temp = amp[int_i_stride + stride] + amp_address[int_j_stride]; amp1 = temp + amp[int_i_stride]; - triangle_interpolation2D(&freq1[i], &freq1[i + 1], &freq1_address[j], &freq2[i], - &freq2[i + 1], &freq2_address[j], &1, spec, m0, m1, + triangle_interpolation2D(freq1[i], freq1[i + 1], freq1_address[j], freq2[i], + freq2[i + 1], freq2_address[j], amp1, spec, m0, m1, iso_intrp); if (i < local_index) { temp += amp_address[int_j_stride + stride]; - triangle_interpolation2D(&freq1[i + 1], &freq1_address[j], &freq1_address[j + 1], - &freq2[i + 1], &freq2_address[j], &freq2_address[j + 1], - &temp, spec, m0, m1, iso_intrp); + triangle_interpolation2D(freq1[i + 1], freq1_address[j], freq1_address[j + 1], + freq2[i + 1], freq2_address[j], freq2_address[j + 1], + temp, spec, m0, m1, iso_intrp); } else { local_index = j + nt; i++; diff --git a/src/c_lib/lib/method.c b/src/c_lib/lib/method.c index 341b1c784..f3cd45047 100644 --- a/src/c_lib/lib/method.c +++ b/src/c_lib/lib/method.c @@ -21,8 +21,8 @@ void MRS_free_event(MRS_event *the_event) { } /** Free the memory/buffer allocation for the MRS dimensions and events within. **/ -void MRS_free_dimension(MRS_dimension *dimensions, unsigned int n) { - unsigned int dim, evt; +void MRS_free_dimension(MRS_dimension *dimensions, int n) { + int dim, evt; for (dim = 0; dim < n; dim++) { if (DEBUG) printf("post execution dimension %d cleanup\n", dim); for (evt = 0; evt < dimensions[dim].n_events; evt++) { @@ -143,31 +143,30 @@ static inline void create_plans_for_events_in_dimension( double coordinates_offset, int n_events, double *fraction, double *duration, unsigned char *is_spectral, double *rotor_frequency_in_Hz, double *rotor_angle_in_rad, double *magnetic_flux_density_in_T, - unsigned int number_of_sidebands) { + int number_of_sidebands) { int i; dim->count = count; dim->coordinates_offset = coordinates_offset; dim->increment = increment; dim->n_events = n_events; - dim->events = (MRS_event *)malloc(n_events * sizeof(MRS_event)); + dim->events = (MRS_event *)malloc((size_t)n_events * sizeof(MRS_event)); dim->inverse_increment = 1.0 / increment; dim->normalize_offset = 0.5 - (coordinates_offset * dim->inverse_increment); dim->R0_offset = 0.0; - MRS_plan *plan = - MRS_create_plan(scheme, number_of_sidebands, *rotor_frequency_in_Hz, - *rotor_angle_in_rad, increment, scheme->allow_4th_rank); + MRS_plan *plan = MRS_create_plan(scheme, number_of_sidebands, *rotor_frequency_in_Hz, + *rotor_angle_in_rad, scheme->allow_4th_rank); // normalize the sideband frequencies. cblas_dscal(number_of_sidebands, dim->inverse_increment, plan->vr_freq, 1); for (i = 0; i < n_events; i++) { - dim->events[i].event_freq_amplitude = malloc_complex128(plan->size); + dim->events[i].event_freq_amplitude = malloc_complex128((size_t)plan->size); vm_double_ones(plan->size * 2, (double *)dim->events[i].event_freq_amplitude); cblas_dscal(plan->size, 0.0, (double *)dim->events[i].event_freq_amplitude + 1, 2); // if (*rotor_frequency_in_Hz != 0.0 && *rotor_frequency_in_Hz != 1.0e12) { - // dim->events[i].freq_amplitude = malloc_double(plan->size); + // dim->events[i].freq_amplitude = malloc_double((size_t)plan->size); // vm_double_ones(plan->size, dim->events[i].freq_amplitude); // } MRS_set_event(&(dim->events[i]), *fraction++, *duration++, *is_spectral++, @@ -181,10 +180,12 @@ static inline void create_plans_for_events_in_dimension( /* buffer to hold the local frequencies and frequency offset. The buffer is useful * when the rotor angle is off magic angle (54.735 deg). */ - dim->local_frequency = malloc_double(scheme->n_gamma * scheme->total_orientations); - dim->local_phase = malloc_double(scheme->n_gamma * scheme->total_orientations); - dim->freq_offset = malloc_double(scheme->octant_orientations); - dim->freq_amplitude = malloc_double(plan->size); + dim->local_frequency = + malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + dim->local_phase = + malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + dim->freq_offset = malloc_double((size_t)scheme->octant_orientations); + dim->freq_amplitude = malloc_double((size_t)plan->size); } /** @@ -195,10 +196,10 @@ MRS_dimension *MRS_create_dimensions( MRS_averaging_scheme *scheme, int *count, double *coordinates_offset, double *increment, double *fractions, double *durations, unsigned char *is_spectral, double *magnetic_flux_density_in_T, double *rotor_frequency_in_Hz, - double *rotor_angle_in_rad, int *n_events, unsigned int n_dim, - unsigned int *number_of_sidebands) { - unsigned int i; - MRS_dimension *dimension = (MRS_dimension *)malloc(n_dim * sizeof(MRS_dimension)); + double *rotor_angle_in_rad, int *n_events, int n_dim, int *number_of_sidebands) { + int i; + MRS_dimension *dimension = + (MRS_dimension *)malloc((size_t)n_dim * sizeof(MRS_dimension)); for (i = 0; i < n_dim; i++) { create_plans_for_events_in_dimension( @@ -213,21 +214,23 @@ MRS_dimension *MRS_create_dimensions( rotor_angle_in_rad += n_events[i]; magnetic_flux_density_in_T += n_events[i]; - // printf("dimension %d\n", i); - // printf("\tcount %d\n", dimension[i].count); - // printf("\tincrement %f Hz\n", dimension[i].increment); - // printf("\tcoordinates offset %f Hz\n", dimension[i].coordinates_offset); - // for (int j = 0; j < n_events[i]; j++) { - // printf("\tEvent %d\n", j); - // printf("\t\tfraction %f\n", dimension[i].events[j].fraction); - // printf("\t\trotor frequency %f Hz\n", - // dimension[i].events[j].rotor_frequency_in_Hz); - // printf("\t\trotor angle %f rad\n", - // dimension[i].events[j].rotor_angle_in_rad); - // printf("\t\tmagnetic flux density %f T\n", - // dimension[i].events[j].magnetic_flux_density_in_T); - // printf("\n"); - // } + if (DEBUG) { + printf("Creating dimensions"); + printf("dimension %d\n", i); + printf("\tcount %d\n", dimension[i].count); + printf("\tincrement %f Hz\n", dimension[i].increment); + printf("\tcoordinates offset %f Hz\n", dimension[i].coordinates_offset); + for (int j = 0; j < n_events[i]; j++) { + printf("\tEvent %d\n", j); + printf("\t\tfraction %f\n", dimension[i].events[j].fraction); + printf("\t\trotor frequency %f Hz\n", + dimension[i].events[j].rotor_frequency_in_Hz); + printf("\t\trotor angle %f rad\n", dimension[i].events[j].rotor_angle_in_rad); + printf("\t\tmagnetic flux density %f T\n", + dimension[i].events[j].magnetic_flux_density_in_T); + printf("\n"); + } + } } return dimension; } diff --git a/src/c_lib/lib/mrsimulator.c b/src/c_lib/lib/mrsimulator.c index f05f48726..8972759b1 100644 --- a/src/c_lib/lib/mrsimulator.c +++ b/src/c_lib/lib/mrsimulator.c @@ -81,10 +81,9 @@ void MRS_plan_release_temp_storage(MRS_plan *the_plan) { * 4) creating the fftw plan, 4) allocating buffer for storing the evaluated frequencies * and their respective amplitudes. */ -MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, - unsigned int number_of_sidebands, +MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, - double increment, bool allow_4th_rank) { + bool allow_4th_rank) { MRS_plan *plan = malloc(sizeof(MRS_plan)); plan->number_of_sidebands = number_of_sidebands; plan->rotor_frequency_in_Hz = rotor_frequency_in_Hz; @@ -113,7 +112,7 @@ MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, * Normalizing amplitudes from the spherical averaging scheme by the number of * sidebands square times the number of octants. */ - plan->norm_amplitudes = malloc_double(scheme->octant_orientations); + plan->norm_amplitudes = malloc_double((size_t)scheme->octant_orientations); cblas_dcopy(scheme->octant_orientations, scheme->amplitudes, 1, plan->norm_amplitudes, 1); double scale = (1.0 / (double)(plan->number_of_sidebands * plan->number_of_sidebands * @@ -135,8 +134,7 @@ MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, */ void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, double rotor_frequency_in_Hz) { - unsigned int size_4; - // double increment_inverse = 1.0 / increment; + int size_4; plan->rotor_frequency_in_Hz = rotor_frequency_in_Hz; plan->is_static = (rotor_frequency_in_Hz < 1.0e-3) ? true : false; @@ -151,7 +149,7 @@ void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, * @see get_sideband_phase_components() */ size_4 = 4 * plan->number_of_sidebands; - plan->pre_phase = malloc_complex128(size_4); + plan->pre_phase = malloc_complex128((size_t)size_4); get_sideband_phase_components(plan->number_of_sidebands, rotor_frequency_in_Hz, (double *)plan->pre_phase); } @@ -165,14 +163,14 @@ void MRS_plan_update_from_rotor_frequency_in_Hz(MRS_plan *plan, */ void MRS_plan_update_from_rotor_angle_in_rad(MRS_plan *plan, double rotor_angle_in_rad, bool allow_4th_rank) { - unsigned int size_2, size_4, i, j; + int size_2, size_4, i, j; plan->rotor_angle_in_rad = rotor_angle_in_rad; /** * Calculate wigner-2j d^2_{m,0} vector where m ∈ [-2, 2]. This vector is used to * rotate the second-rank tensors from the rotor frame to the lab frame. * @see wigner_dm0_vector() */ - plan->wigner_d2m0_vector = malloc_double(5); + plan->wigner_d2m0_vector = malloc_double((size_t)5); wigner_dm0_vector(2, rotor_angle_in_rad, plan->wigner_d2m0_vector); plan->wigner_d4m0_vector = NULL; @@ -182,13 +180,13 @@ void MRS_plan_update_from_rotor_angle_in_rad(MRS_plan *plan, double rotor_angle_ * rotate the fourth-rank tensors from the rotor frame to the lab frame. * @see wigner_dm0_vector() */ - plan->wigner_d4m0_vector = malloc_double(9); + plan->wigner_d4m0_vector = malloc_double((size_t)9); wigner_dm0_vector(4, rotor_angle_in_rad, plan->wigner_d4m0_vector); } // pre_phase_2 is only calculated for m=-2 and -1 for l=2 rank tensor calculation. size_2 = 2 * plan->number_of_sidebands; - plan->pre_phase_2 = malloc_complex128(size_2); + plan->pre_phase_2 = malloc_complex128((size_t)size_2); /* Copy the pre_phase[m=-2 to 2] to pre_phase2 */ cblas_zcopy(size_2, (double *)(plan->pre_phase[2 * plan->number_of_sidebands]), 1, @@ -216,7 +214,7 @@ void MRS_plan_update_from_rotor_angle_in_rad(MRS_plan *plan, double rotor_angle_ /* pre_phase_4 is only calculated for m=-4, -3, -2, and -1 for l=4 rank tensor * calculation. */ size_4 = 4 * plan->number_of_sidebands; - plan->pre_phase_4 = malloc_complex128(size_4); + plan->pre_phase_4 = malloc_complex128((size_t)size_4); /* Copy the pre_phase[m=-4 to 4] to pre_phase4 */ cblas_zcopy(size_4, (double *)(plan->pre_phase), 1, (double *)(plan->pre_phase_4), 1); @@ -406,8 +404,8 @@ void MRS_get_normalized_frequencies_from_plan(MRS_averaging_scheme *scheme, double fraction, unsigned char is_spectral, double duration) { - unsigned int i, g_idx, ptr; - unsigned int freq_size = scheme->n_gamma * scheme->total_orientations; + int i, g_idx, ptr; + int freq_size = scheme->n_gamma * scheme->total_orientations; double temp, fraction_duration; double *f_complex, *local_frequency; @@ -442,7 +440,7 @@ void MRS_get_normalized_frequencies_from_plan(MRS_averaging_scheme *scheme, /* Normalized local anisotropic frequency contributions from the 2nd-rank tensor. */ for (g_idx = 0; g_idx < scheme->n_gamma; g_idx++) { ptr = scheme->total_orientations * g_idx; - temp = 2 * fraction_duration; + temp = 2.0 * fraction_duration; if (plan->is_static) { for (i = 0; i < 2; i++) { @@ -508,7 +506,7 @@ static inline void MRS_rotate_single_site_interaction_components( double B0_in_T, // Magnetic flux density in T. unsigned char *freq_contrib // The pointer to freq contribs boolean. ) { - unsigned int i, n_sites = sites->number_of_sites; + int i, n_sites = sites->number_of_sites; double larmor_freq_in_MHz, larmor_freq_in_Hz; float *mf = &transition[n_sites], *mi = transition; double F2_shield[10]; @@ -588,14 +586,14 @@ static inline void quad_coupling_cross_terms( coupling_struct *couplings, // Pointer to a list of couplings within a spin system. site_struct *sites, // Pointer to a list of sites within a spin system. int site_index, // site index - unsigned int coupling_index, // coupled index - double *F0, // The F0 components. - complex128 *F2, // The F2 components. - complex128 *F4, // The F4 components. - double *F0_temp, // The temporary F0 components. - complex128 *F2_temp, // The temporary F2 components. - complex128 *F4_temp, // The temporary F4 components. - double B0_in_T, // Magnetic flux density in T. + int coupling_index, // coupled index + double *F0, // The F0 components. + complex128 *F2, // The F2 components. + complex128 *F4, // The F4 components. + double *F0_temp, // The temporary F0 components. + complex128 *F2_temp, // The temporary F2 components. + complex128 *F4_temp, // The temporary F4 components. + double B0_in_T, // Magnetic flux density in T. unsigned char *freq_contrib, // The pointer to freq contribs boolean. float mIi, // inital transition state of site coupled to the quad site (p) float mSqi, // inital transition state of the quad site (d) @@ -642,7 +640,7 @@ static inline void MRS_rotate_coupled_site_interaction_components( coupling_struct *couplings, // Pointer to a list of couplings within a spin system. site_struct *sites, // Pointer to a list of sites within a spin system. float *transition, // The spin transition. - unsigned int n_sites, // The number of sites. + int n_sites, // The number of sites. double *F0, // The F0 components. complex128 *F2, // The F2 components. complex128 *F4, // The F4 components. @@ -652,7 +650,7 @@ static inline void MRS_rotate_coupled_site_interaction_components( double B0_in_T, // Magnetic flux density in T. unsigned char *freq_contrib // The pointer to freq contribs boolean. ) { - unsigned int i, j = 0, n_couplings = couplings->number_of_couplings; + int i, j = 0, n_couplings = couplings->number_of_couplings; int site_index_I, site_index_S; float mIf, mSf, mIi, mSi; @@ -759,15 +757,15 @@ void MRS_rotate_components_from_PAS_to_common_frame( * = scale [-[cos(m ωr t) -1] +Isin(m ωr t)] * That is, pre_phase[-m] = -Re(pre_phase[m]) + Im(pre_phase[m]) */ -void get_sideband_phase_components_2(unsigned int number_of_sidebands, +void get_sideband_phase_components_2(int number_of_sidebands, double rotor_frequency_in_Hz, complex128 *pre_phase) { int m, i; double spin_angular_freq, tau, scale; - double *input = malloc_double(number_of_sidebands); - double *ones = malloc_double(number_of_sidebands); - double *phase = malloc_double(number_of_sidebands); + double *input = malloc_double((size_t)number_of_sidebands); + double *ones = malloc_double((size_t)number_of_sidebands); + double *phase = malloc_double((size_t)number_of_sidebands); vm_double_ones(number_of_sidebands, ones); vm_double_arange(number_of_sidebands, input); @@ -835,11 +833,11 @@ void get_sideband_phase_components_2(unsigned int number_of_sidebands, * as the leading dimension. The first number_of_sidebands entries corresponds to * m_wr=-4. */ -void get_sideband_phase_components(unsigned int number_of_sidebands, +void get_sideband_phase_components(int number_of_sidebands, double sample_rotation_frequency, double *restrict pre_phase) { double spin_angular_freq, tau, wrt, pht, scale; - unsigned int step, m; + int step, m; // Calculate the spin angular frequency spin_angular_freq = sample_rotation_frequency * CONST_2PI; diff --git a/src/c_lib/lib/octahedron.c b/src/c_lib/lib/octahedron.c index 3ebe0f834..72dc6e746 100644 --- a/src/c_lib/lib/octahedron.c +++ b/src/c_lib/lib/octahedron.c @@ -9,12 +9,12 @@ #include "octahedron.h" -void octahedronGetDirectionCosineSquareAndWeightsOverOctant(const unsigned int nt, +void octahedronGetDirectionCosineSquareAndWeightsOverOctant(const int nt, double *restrict xr, double *restrict yr, double *restrict zr, double *restrict amp) { - unsigned int i, j; + int i, j; // scale is divided by pi to remove angular dependence. // The factor 6 is because each point in the triangle interpolation is shared by // 6 neighboring triangles. @@ -48,15 +48,14 @@ void octahedronGetDirectionCosineSquareAndWeightsOverOctant(const unsigned int n *amp = 1.0 / (r2 * r2 * factor); } -void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, - double *restrict cos_alpha, +void octahedronGetPolarAngleTrigOverOctant(const int nt, double *restrict cos_alpha, double *restrict cos_beta, double *restrict amp) { - unsigned int points = (nt + 1) * (nt + 2) / 2; - double *xr = malloc_double(points); - double *yr = malloc_double(points); - double *zr = malloc_double(points); - double *sin_beta = malloc_double(points); + int points = (nt + 1) * (nt + 2) / 2; + double *xr = malloc_double((size_t)points); + double *yr = malloc_double((size_t)points); + double *zr = malloc_double((size_t)points); + double *sin_beta = malloc_double((size_t)points); // The values xr = x^2, yr = y^2, zr = z^2, where x, y, and z are the // direction cosines. @@ -78,7 +77,7 @@ void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, // vdSqrt(points, &yr[0], &yr[0]); // Evaluate cos_alpha = x/sqrt(x^2 + y^2) = zr/sin_beta - vm_double_divide(points - 1, zr, sin_beta, cos_alpha); + vm_double_divide((points - 1), zr, sin_beta, cos_alpha); // vdDiv(points-1, yr, sinBeta, sinAlpha ); cos_alpha[points - 1] = 1.0; @@ -90,14 +89,14 @@ void octahedronGetPolarAngleTrigOverOctant(const unsigned int nt, free(sin_beta); } -void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, +void octahedronGetComplexExpOfPolarAngleOverOctant(const int nt, void *restrict exp_I_alpha, void *restrict exp_I_beta, double *restrict amp) { - unsigned int points = (nt + 1) * (nt + 2) / 2; - double *xr = malloc_double(points); - double *yr = malloc_double(points); - double *zr = malloc_double(points); + int points = (nt + 1) * (nt + 2) / 2; + double *xr = malloc_double((size_t)points); + double *yr = malloc_double((size_t)points); + double *zr = malloc_double((size_t)points); octahedronGetDirectionCosineSquareAndWeightsOverOctant(nt, xr, yr, zr, amp); // At this point the variables @@ -143,7 +142,7 @@ void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, // cos(alpha) = x/sqrt(x^2 + y^2) = x/sin(beta) ... (3) // In terms of the variables, Eq (3) is given as xr/zr. // Evaluate xr = xr/zr - vm_double_divide_inplace(points - 1, zr, xr); + vm_double_divide_inplace((points - 1), zr, xr); xr[points - 1] = 0.0; // Copy cos(alpha), aka xr, to the even addresses of exp_I_alpha cblas_dcopy(points, xr, 1, (double *)exp_I_alpha, 2); @@ -152,7 +151,7 @@ void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, // sin(alpha) = y/sqrt(x^2 + y^2) = y/sin(beta) ... (4) // In terms of the variables, Eq (4) is given as yr/zr. // Evaluate yr = yr/zr - vm_double_divide_inplace(points - 1, zr, yr); + vm_double_divide_inplace((points - 1), zr, yr); yr[points - 1] = 0.0; // Copy sin(alpha), aka yr, to the odd addresses of exp_I_alpha. cblas_dcopy(points, yr, 1, (double *)exp_I_alpha + 1, 2); @@ -163,8 +162,8 @@ void octahedronGetComplexExpOfPolarAngleOverOctant(const unsigned int nt, free(zr); } -void get_total_amplitude(const unsigned int nt, double *amp, double *amp_sum) { - unsigned int i = 0, j = 0, local_index = nt - 1, n_pts = (nt + 1) * (nt + 2) / 2; +void get_total_amplitude(const int nt, double *amp, double *amp_sum) { + int i = 0, j = 0, local_index = nt - 1, n_pts = (nt + 1) * (nt + 2) / 2; double *amp_address, temp; amp_address = &[nt + 1]; @@ -189,8 +188,8 @@ void get_total_amplitude(const unsigned int nt, double *amp, double *amp_sum) { // fix amplitude for binning case: // octant face edge points are divided by 2 // octant vertexes are divided by 4 -void fix_amplitude_for_binning(unsigned int nt, double *amp) { - unsigned int i = 0, i_max = (nt + 1) * (nt + 2) / 2, tt = nt - 1; +void fix_amplitude_for_binning(int nt, double *amp) { + int i = 0, i_max = (nt + 1) * (nt + 2) / 2, tt = nt - 1; bool inc_one = false; while (i <= nt) amp[i++] /= 2.0; while (i < i_max) { @@ -211,7 +210,7 @@ void fix_amplitude_for_binning(unsigned int nt, double *amp) { cblas_dscal(i_max, 6.0, amp, 1); } -void averaging_setup(unsigned int nt, void *exp_I_alpha, void *exp_I_beta, double *amp, +void averaging_setup(int nt, void *exp_I_alpha, void *exp_I_beta, double *amp, bool interpolation) { // octahedronGetPolarAngleTrigOverOctant(nt, cos_alpha, cos_beta, amp); octahedronGetComplexExpOfPolarAngleOverOctant(nt, exp_I_alpha, exp_I_beta, amp); diff --git a/src/c_lib/lib/schemes.c b/src/c_lib/lib/schemes.c index 50c343fa3..3f8485661 100644 --- a/src/c_lib/lib/schemes.c +++ b/src/c_lib/lib/schemes.c @@ -11,7 +11,7 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, complex128 *exp_I_beta, bool allow_4th_rank) { - unsigned int allocate_size_2, allocate_size_4; + int allocate_size_2, allocate_size_4; double delta_alpha = 0.0; scheme->total_orientations = scheme->octant_orientations; @@ -56,7 +56,7 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, /* Second-rank reduced wigner matrices at every β orientation from the positive upper * octant. */ - scheme->wigner_2j_matrices = malloc_double(allocate_size_2); + scheme->wigner_2j_matrices = malloc_double((size_t)allocate_size_2); wigner_d_matrices_from_exp_I_beta(2, scheme->octant_orientations, true, exp_I_beta, scheme->wigner_2j_matrices); @@ -64,7 +64,7 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, if (allow_4th_rank) { /* Fourt-rank reduced wigner matrices at every β orientation from the positive upper * octant. */ - scheme->wigner_4j_matrices = malloc_double(allocate_size_4); + scheme->wigner_4j_matrices = malloc_double((size_t)allocate_size_4); wigner_d_matrices_from_exp_I_beta(4, scheme->octant_orientations, true, exp_I_beta, scheme->wigner_4j_matrices); } @@ -106,13 +106,13 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, /* ................................................................................ */ /* w2 is the buffer for storing the frequencies calculated from the second-rank * tensors. Only calcuate the -2, -1, and 0 tensor components.*/ - scheme->w2 = malloc_complex128(3 * scheme->total_orientations); + scheme->w2 = malloc_complex128((size_t)(3 * scheme->total_orientations)); scheme->w4 = NULL; if (allow_4th_rank) { /* w4 is the buffer for storing the frequencies calculated from the fourth-rank * tensors. Only calcuate the -4, -3, -2, -1, and 0 tensor components.*/ - scheme->w4 = malloc_complex128(5 * scheme->total_orientations); + scheme->w4 = malloc_complex128((size_t)(5 * scheme->total_orientations)); } } @@ -135,10 +135,9 @@ void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) { } /* Create a new orientation averaging scheme. */ -MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_density, - bool allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, +MRS_averaging_scheme *MRS_create_averaging_scheme(int integration_density, + bool allow_4th_rank, int n_gamma, + int integration_volume, bool interpolation) { int alpha_size; MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); @@ -160,10 +159,10 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi // The 4 * octant_orientations memory allocation is for m=4, 3, 2, and 1 alpha_size = 4 * scheme->octant_orientations; alpha_size *= (integration_volume == 2) ? 2 : 1; - scheme->exp_Im_alpha = malloc_complex128(alpha_size); - scheme->exp_Im_gamma = malloc_complex128(4 * n_gamma); - complex128 *exp_I_beta = malloc_complex128(scheme->octant_orientations); - scheme->amplitudes = malloc_double(scheme->octant_orientations); + scheme->exp_Im_alpha = malloc_complex128((size_t)alpha_size); + scheme->exp_Im_gamma = malloc_complex128((size_t)(4 * n_gamma)); + complex128 *exp_I_beta = malloc_complex128((size_t)scheme->octant_orientations); + scheme->amplitudes = malloc_double((size_t)scheme->octant_orientations); averaging_setup(integration_density, &scheme->exp_Im_alpha[3 * scheme->octant_orientations], exp_I_beta, @@ -172,8 +171,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi averaging_scheme_setup(scheme, exp_I_beta, allow_4th_rank); // exp(-im gamma) for m=[-4,-1], and gamma=[0..8]*2pi/9 - double *gamma = malloc_double(n_gamma); - double *temp = malloc_double(n_gamma); + double *gamma = malloc_double((size_t)n_gamma); + double *temp = malloc_double((size_t)n_gamma); double factor = CONST_2PI / (double)n_gamma; vm_double_arange(n_gamma, gamma); cblas_dscal(n_gamma, factor, gamma, 1); @@ -182,10 +181,11 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi // reallocate exp_I_beta memory as scratch. scheme->scratch = (double *)exp_I_beta; - scheme->amps_real = malloc_double(scheme->total_orientations); - scheme->amps_imag = malloc_double(scheme->total_orientations); - scheme->phase = malloc_double(scheme->n_gamma * scheme->total_orientations); - scheme->exp_I_phase = malloc_complex128(scheme->n_gamma * scheme->total_orientations); + scheme->amps_real = malloc_double((size_t)scheme->total_orientations); + scheme->amps_imag = malloc_double((size_t)scheme->total_orientations); + scheme->phase = malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + scheme->exp_I_phase = + malloc_complex128((size_t)(scheme->n_gamma * scheme->total_orientations)); free(gamma); free(temp); @@ -194,9 +194,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme(unsigned int integration_densi /* Create a new orientation averaging scheme. */ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( - double *alpha, double *beta, double *weight, unsigned int n_angles, - bool allow_4th_rank, unsigned int n_gamma, const unsigned int position_size, - int32_t *positions, bool interpolation) { + double *alpha, double *beta, double *weight, int n_angles, bool allow_4th_rank, + int n_gamma, const int position_size, int32_t *positions, bool interpolation) { double scale; MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); @@ -214,10 +213,10 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( /* Calculate α, β, and weights over the positive octant. .......................... */ /* ................................................................................ */ // The 4 * octant_orientations memory allocation is for m=4, 3, 2, and 1 - scheme->exp_Im_alpha = malloc_complex128(4 * scheme->octant_orientations); - scheme->exp_Im_gamma = malloc_complex128(4 * n_gamma); - complex128 *exp_I_beta = malloc_complex128(scheme->octant_orientations); - scheme->amplitudes = malloc_double(scheme->octant_orientations); + scheme->exp_Im_alpha = malloc_complex128((size_t)(4 * scheme->octant_orientations)); + scheme->exp_Im_gamma = malloc_complex128((size_t)(4 * n_gamma)); + complex128 *exp_I_beta = malloc_complex128((size_t)scheme->octant_orientations); + scheme->amplitudes = malloc_double((size_t)scheme->octant_orientations); // scale = (1 / 6) factor for 1 point to 6 triangle ratio. scale = (interpolation) ? 0.1666666667 : 1.0; @@ -233,8 +232,8 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( averaging_scheme_setup(scheme, exp_I_beta, allow_4th_rank); // exp(-im gamma) for m=[-4,-1], and gamma=[0..8]*2pi/9 - double *gamma = malloc_double(n_gamma); - double *temp = malloc_double(n_gamma); + double *gamma = malloc_double((size_t)n_gamma); + double *temp = malloc_double((size_t)n_gamma); double factor = CONST_2PI / (double)n_gamma; vm_double_arange(n_gamma, gamma); cblas_dscal(n_gamma, factor, gamma, 1); @@ -243,10 +242,11 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( // reallocate exp_I_beta memory as scratch. scheme->scratch = (double *)exp_I_beta; - scheme->amps_real = malloc_double(scheme->total_orientations); - scheme->amps_imag = malloc_double(scheme->total_orientations); - scheme->phase = malloc_double(scheme->n_gamma * scheme->total_orientations); - scheme->exp_I_phase = malloc_complex128(scheme->n_gamma * scheme->total_orientations); + scheme->amps_real = malloc_double((size_t)scheme->total_orientations); + scheme->amps_imag = malloc_double((size_t)scheme->total_orientations); + scheme->phase = malloc_double((size_t)(scheme->n_gamma * scheme->total_orientations)); + scheme->exp_I_phase = + malloc_complex128((size_t)(scheme->n_gamma * scheme->total_orientations)); free(gamma); free(temp); @@ -256,15 +256,15 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( /* ---------------------------------------------------------------------------------- */ /* fftw routine setup ............................................................... */ /* .................................................................................. */ -MRS_fftw_scheme *create_fftw_scheme(unsigned int total_orientations, - unsigned int number_of_sidebands) { - unsigned int size = total_orientations * number_of_sidebands; - int nssb = (int)number_of_sidebands; +MRS_fftw_scheme *create_fftw_scheme(int total_orientations, int number_of_sidebands) { + int size = total_orientations * number_of_sidebands; + int nssb = number_of_sidebands; MRS_fftw_scheme *fftw_scheme = malloc(sizeof(MRS_fftw_scheme)); // fftw_scheme->vector = fftw_alloc_complex(size); - fftw_scheme->vector = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * size); - // malloc_complex128(plan->size); + fftw_scheme->vector = + (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (size_t)size); + // malloc_complex128((size_t)plan->size); // gettimeofday(&fft_setup_time, NULL); // int fftw_thread = fftw_init_threads(); diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index cf3a5b699..a36096087 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -50,7 +50,7 @@ void __mrsimulator_core( MRS_dimension *dimensions, // Pointer to MRS_dimension structure. MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. - unsigned int iso_intrp, // Isotropic interpolation scheme (linear | Gaussian) + int iso_intrp, // Isotropic interpolation scheme (linear | Gaussian) unsigned char *freq_contrib, // A list of freq_contrib booleans. double *affine_matrix // Affine transformation matrix. ) { @@ -59,8 +59,9 @@ void __mrsimulator_core( et al. `Computation of Orientational Averages in Solid-State NMR by Gaussian Spherical Quadrature` JMR, 132, 1998. https://doi.org/10.1006/jmre.1998.1427 */ + unsigned char is_spectral; - unsigned int evt; + int evt; int dim, total_pts = scheme->n_gamma * scheme->total_orientations; double B0_in_T, fraction, duration; vm_double_zeros(total_pts, scheme->phase); @@ -156,10 +157,10 @@ void __mrsimulator_core( void mrsimulator_core( // spectrum information and related amplitude - double *spec, // Pointer to the spectrum array (complex). - double coordinates_offset, // The start of the frequency spectrum. - double increment, // The increment of the frequency spectrum. - int count, // Number of points on the frequency spectrum. + double *spec, // Pointer to the spectrum array (complex). + // double coordinates_offset, // The start of the frequency spectrum. + // double increment, // The increment of the frequency spectrum. + // int count, // Number of points on the frequency spectrum. site_struct *sites, // Pointer to a list of sites within a spin system. coupling_struct *couplings, // Pointer to a list of couplings within a spin system. MRS_dimension *dimensions, // the dimensions in the method. @@ -167,16 +168,16 @@ void mrsimulator_core( int quad_second_order, // Quad theory for second order, // spin rate, spin angle and number spinning sidebands - unsigned int number_of_sidebands, // The number of sidebands - double rotor_frequency_in_Hz, // The rotor spin frequency - double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis + int number_of_sidebands, // The number of sidebands + double rotor_frequency_in_Hz, // The rotor spin frequency + double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis // Pointer to the a list of transitions. float *transition_pathway, double *transition_pathway_weight, // powder orientation average int integration_density, // The number of triangle along the edge of octahedron - unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. - unsigned int interpolate_type, unsigned char *freq_contrib, double *affine_matrix) { + int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. + int interpolate_type, unsigned char *freq_contrib, double *affine_matrix) { // int num_process = openblas_get_num_procs(); // int num_threads = openblas_get_num_threads(); // openblas_set_num_threads(1); diff --git a/src/c_lib/sandbox/sandbox.pxd b/src/c_lib/sandbox/sandbox.pxd index ff2d672de..cd67c92fb 100644 --- a/src/c_lib/sandbox/sandbox.pxd +++ b/src/c_lib/sandbox/sandbox.pxd @@ -12,22 +12,22 @@ from libc.stdint cimport int32_t cdef extern from "schemes.h": ctypedef struct MRS_averaging_scheme: - unsigned int total_orientations - unsigned int integration_density - unsigned int integration_volume + int total_orientations + int integration_density + int integration_volume MRS_averaging_scheme * MRS_create_averaging_scheme( - unsigned int integration_density, + int integration_density, bool_t allow_4th_rank, - unsigned int n_gamma, - unsigned int integration_volume, + int n_gamma, + int integration_volume, bool_t interpolation) MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( double *alpha, double *beta, - double *weight, unsigned int n_angles, + double *weight, int n_angles, bool_t allow_4th_rank, - const unsigned int position_size, + const int position_size, int32_t *positions, bool_t interpolation) @@ -36,13 +36,13 @@ cdef extern from "schemes.h": cdef extern from "mrsimulator.h": ctypedef struct MRS_plan: - unsigned int number_of_sidebands + int number_of_sidebands double rotor_frequency_in_Hz double rotor_angle_in_rad - MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, unsigned int number_of_sidebands, + MRS_plan *MRS_create_plan(MRS_averaging_scheme *scheme, int number_of_sidebands, double rotor_frequency_in_Hz, - double rotor_angle_in_rad, double increment, + double rotor_angle_in_rad, bool_t allow_4th_rank) void MRS_free_plan(MRS_plan *plan) void MRS_get_amplitudes_from_plan(MRS_plan *plan, double complex *R2, diff --git a/src/c_lib/sandbox/sandbox.pyx b/src/c_lib/sandbox/sandbox.pyx index cf0d9efce..519c3b0d1 100644 --- a/src/c_lib/sandbox/sandbox.pyx +++ b/src/c_lib/sandbox/sandbox.pyx @@ -168,7 +168,7 @@ cdef class MRSPlan: self.plan = clib.MRS_create_plan(averaging_scheme.scheme, number_of_sidebands, rotor_frequency_in_Hz, - rotor_angle_in_rad, increment, allow_4th_rank) + rotor_angle_in_rad, allow_4th_rank) @property def number_of_sidebands(self): diff --git a/src/c_lib/test/test.pxd b/src/c_lib/test/test.pxd index 79acd738f..ad6b7be44 100644 --- a/src/c_lib/test/test.pxd +++ b/src/c_lib/test/test.pxd @@ -65,11 +65,11 @@ cdef extern from "angular_momentum/wigner_matrix.h": void wigner_dm0_vector(const int l, const double beta, double *R_out) - void get_exp_Im_angle(const unsigned int octant_orientations, + void get_exp_Im_angle(const int octant_orientations, const bool_t allow_4th_rank, void *exp_Im_angle, double delta_alpha) - void __batch_wigner_rotation(const unsigned int octant_orientations, - const unsigned int n_octants, double *wigner_2j_matrices, void *R2, + void __batch_wigner_rotation(const int octant_orientations, + const int n_octants, double *wigner_2j_matrices, void *R2, double *wigner_4j_matrices, void *R4, void *exp_Im_alpha, void *w2, void *w4) @@ -84,64 +84,64 @@ cdef extern from "octahedron.h": cdef extern from "interpolation.h": void triangle_interpolation1D( - double *freq1, - double *freq2, - double *freq3, - double *amp, + double freq1, + double freq2, + double freq3, + double amp, double *spec, - int *points, - unsigned int iso_intrp) + int points, + int iso_intrp) void triangle_interpolation1D_linear( - double *freq1, - double *freq2, - double *freq3, - double *amp, + double freq1, + double freq2, + double freq3, + double amp, double *spec, - int *points) + int points) void triangle_interpolation1D_gaussian( - double *freq1, - double *freq2, - double *freq3, - double *amp, + double freq1, + double freq2, + double freq3, + double amp, double *spec, - int *points) + int points) void triangle_interpolation2D( - double *freq11, - double *freq12, - double *freq13, - double *freq21, - double *freq22, - double *freq23, - double *amp, + double freq11, + double freq12, + double freq13, + double freq21, + double freq22, + double freq23, + double amp, double *spec, int m0, int m1, - unsigned int iso_intrp) + int iso_intrp) void octahedronInterpolation( double *spec, double *freq, int nt, double *amp, - int stride, + const int stride, int m) cdef extern from "mrsimulator.h": void get_sideband_phase_components( - unsigned int number_of_sidebands, + int number_of_sidebands, double spin_frequency, double *pre_phase) # ctypedef struct MRS_plan # MRS_plan *MRS_create_plan( -# unsigned int integration_density, +# int integration_density, # int number_of_sidebands, # double rotor_frequency_in_Hz, -# double rotor_angle_in_rad, double increment, +# double rotor_angle_in_rad, # bool_t allow_4th_rank) @@ -173,7 +173,7 @@ cdef extern from "method.h": double increment # Increment of coordinates along the dimension. double coordinates_offset # Start coordinate of the dimension. MRS_event *events # Holds a list of events. - unsigned int n_events # The number of events. + int n_events # The number of events. # MRS_dimension *MRS_create_dimensions( # MRS_averaging_scheme *scheme, @@ -183,7 +183,7 @@ cdef extern from "method.h": # double *magnetic_flux_density_in_T, # double *rotor_frequency_in_Hz, # double *rotor_angle_in_rad, - # unsigned int n_events, + # int n_events, # int number_of_sidebands) cdef extern from "simulation.h": @@ -202,14 +202,14 @@ cdef extern from "simulation.h": # second order quad Hamiltonian. # spin rate, spin angle and number spinning sidebands - unsigned int number_of_sidebands, + int number_of_sidebands, double rotor_frequency_in_Hz, double rotor_angle_in_rad, # The transition as transition[0] = mi and transition[1] = mf double *transition, int integration_density, - unsigned int integration_volume, # 0-octant, 1-hemisphere, 2-sphere. + int integration_volume, # 0-octant, 1-hemisphere, 2-sphere. ) cdef extern from "frequency/spatial_orientation_tensor_components.h": diff --git a/src/c_lib/test/test.pyx b/src/c_lib/test/test.pyx index ba6bc1f11..2d7a3f690 100644 --- a/src/c_lib/test/test.pyx +++ b/src/c_lib/test/test.pyx @@ -23,7 +23,7 @@ def wigner_d_element(float l, float m1, float m2, double beta): @cython.boundscheck(False) @cython.wraparound(False) def wigner_d_matrices(int l, np.ndarray[double] angle): - cdef int n = angle.size + cdef int n = int(angle.size) cdef int n1 = (2*l+1)**2 cdef np.ndarray[double] wigner = np.empty(n*n1, dtype=np.float64) clib.wigner_d_matrices(l, n, &angle[0], &wigner[0]) @@ -46,7 +46,7 @@ def wigner_d_matrices_from_exp_I_beta(int l, bool_t half, np.ndarray[double comp :ivar exp_I_beta: An 1D numpy array or a scalar representing $\exp\beta$. """ n1 = (2 * l + 1) - cdef int n = exp_I_beta.size + cdef int n = int(exp_I_beta.size) size = n * n1*(l+1) if half else n * n1**2 cdef np.ndarray[double, ndim=1] wigner = np.empty(size) @@ -88,7 +88,7 @@ def __wigner_rotation_2(int l, np.ndarray[double] cos_alpha, np.ndarray[double complex] R_in): cdef int n1 = 2 * l + 1 - cdef int n = cos_alpha.size + cdef int n = int(cos_alpha.size) cdef np.ndarray[double, ndim=1] wigner cdef np.ndarray[double complex, ndim=1] exp_I_beta wigner = np.empty(n1 * (l+1) * n, dtype=np.float64) @@ -111,7 +111,7 @@ def __wigner_rotation_2(int l, np.ndarray[double] cos_alpha, @cython.boundscheck(False) @cython.wraparound(False) def get_exp_Im_angle(int n, np.ndarray[double] cos_alpha, bool_t allow_4th_rank): - cdef unsigned int n_ = n + cdef int n_ = n cdef np.ndarray[double complex] exp_Im_angle = np.empty(4*n, dtype=np.complex128) exp_Im_angle[3*n:] = cos_alpha + 1j*np.sqrt(1.0 - cos_alpha**2) clib.get_exp_Im_angle(n_, allow_4th_rank, &exp_Im_angle[0], 0.0) @@ -120,7 +120,7 @@ def get_exp_Im_angle(int n, np.ndarray[double] cos_alpha, bool_t allow_4th_rank) @cython.boundscheck(False) @cython.wraparound(False) -def pre_phase_components(unsigned int number_of_sidebands, double rotor_frequency_in_Hz): +def pre_phase_components(int number_of_sidebands, double rotor_frequency_in_Hz): cdef int n1 = 4 * number_of_sidebands cdef np.ndarray[double] pre_phase = np.zeros(2*n1, dtype=np.float64) clib.get_sideband_phase_components(number_of_sidebands, rotor_frequency_in_Hz, &pre_phase[0]) @@ -157,7 +157,7 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): :return amp: The amplitude at the given $\alpha$ and $\beta$. """ nt = integration_density - cdef unsigned int octant_orientations = int((nt+1) * (nt+2)/2) + cdef int octant_orientations = int((nt+1) * (nt+2)/2) cdef bool_t interploation = True cdef np.ndarray[double complex] exp_I_alpha = np.empty(octant_orientations, dtype=np.complex128) @@ -171,17 +171,17 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): @cython.boundscheck(False) @cython.wraparound(False) -def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, int stride=1): +def octahedronInterpolation(np.ndarray[double] spec, np.ndarray[double, ndim=2] freq, int nt, np.ndarray[double, ndim=2] amp, const int stride=1): cdef int i - cdef int number_of_sidebands = amp.shape[0] + cdef int number_of_sidebands = int(amp.shape[0]) for i in range(number_of_sidebands): - clib.octahedronInterpolation(&spec[0], &freq[i,0], nt, &[i,0], stride, spec.size) + clib.octahedronInterpolation(&spec[0], &freq[i,0], nt, &[i,0], stride, int(spec.size)) @cython.boundscheck(False) @cython.wraparound(False) def triangle_interpolation1D(vector, np.ndarray[double, ndim=1] spectrum_amp, - double amp=1, str type="linear"): + double amp=1.0, str type="linear"): r"""Given a vector of three points, this method interpolates the between the points to form a triangle. The height of the triangle is given as `2.0/(f[2]-f[1])` where `f` is the array `vector` sorted in an ascending @@ -197,24 +197,22 @@ def triangle_interpolation1D(vector, np.ndarray[double, ndim=1] spectrum_amp, default value is 0. :ivar type: Linear or Gaussian interpolation for delta functions. """ - cdef np.ndarray[int, ndim=1] points = np.asarray([spectrum_amp.size/2], dtype=np.int32) + cdef int points = np.asarray([spectrum_amp.size/2], dtype=np.intc)[0] cdef np.ndarray[double, ndim=1] f_vector = np.asarray(vector, dtype=np.float64) - cdef double *f1 = &f_vector[0] - cdef double *f2 = &f_vector[1] - cdef double *f3 = &f_vector[2] - - cdef np.ndarray[double, ndim=1] amp_ = np.asarray([amp]) + cdef double f1 = f_vector[0] + cdef double f2 = f_vector[1] + cdef double f3 = f_vector[2] iso_intrp = 0 if type == "linear" else 1 - clib.triangle_interpolation1D(f1, f2, f3, &_[0], &spectrum_amp[0], - &points[0], iso_intrp) + clib.triangle_interpolation1D(f1, f2, f3, amp, &spectrum_amp[0], + points, iso_intrp) @cython.boundscheck(False) @cython.wraparound(False) def triangle_interpolation2D(vector1, vector2, np.ndarray[double, ndim=2] spectrum_amp, - double amp=1, str type="linear"): + double amp=1.0, str type="linear"): r"""Given a vector of three points, this method interpolates the between the points to form a triangle. The height of the triangle is given as `2.0/(f[2]-f[1])` where `f` is the array `vector` sorted in an ascending @@ -224,29 +222,27 @@ def triangle_interpolation2D(vector1, vector2, np.ndarray[double, ndim=2] spectr :ivar vector2: 1-D array of three points. :ivar spectrum_amp: A numpy array of amplitudes. This array is the output. """ - shape = np.asarray([spectrum_amp.shape[0], spectrum_amp.shape[1]/2], dtype=np.int32) + shape = np.asarray([spectrum_amp.shape[0], spectrum_amp.shape[1]/2], dtype=np.intc) # cdef np.ndarray[int, ndim=1] points = shape cdef np.ndarray[double, ndim=1] f1_vector = np.asarray(vector1, dtype=np.float64) cdef np.ndarray[double, ndim=1] f2_vector = np.asarray(vector2, dtype=np.float64) - cdef double *f11 = &f1_vector[0] - cdef double *f12 = &f1_vector[1] - cdef double *f13 = &f1_vector[2] - - cdef double *f21 = &f2_vector[0] - cdef double *f22 = &f2_vector[1] - cdef double *f23 = &f2_vector[2] + cdef double f11 = f1_vector[0] + cdef double f12 = f1_vector[1] + cdef double f13 = f1_vector[2] - cdef np.ndarray[double, ndim=1] amp_ = np.asarray([amp]) + cdef double f21 = f2_vector[0] + cdef double f22 = f2_vector[1] + cdef double f23 = f2_vector[2] iso_intrp = 0 if type == "linear" else 1 - clib.triangle_interpolation2D(f11, f12, f13, f21, f22, f23, &_[0], + clib.triangle_interpolation2D(f11, f12, f13, f21, f22, f23, amp, &spectrum_amp[0, 0], shape[0], shape[1], iso_intrp) @cython.boundscheck(False) @cython.wraparound(False) -def __batch_wigner_rotation(unsigned int octant_orientations, - unsigned int n_octants, +def __batch_wigner_rotation(int octant_orientations, + int n_octants, np.ndarray[double] wigner_2j_matrices, np.ndarray[double complex] R2, np.ndarray[double] wigner_4j_matrices, @@ -283,122 +279,122 @@ def vm_absd(double a): @cython.wraparound(False) def vm_add(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_add(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_add(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_add_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_add_inplace(A.size, &A[0], &B[0]) + clib.test_vm_double_add_inplace(int(A.size), &A[0], &B[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_sub(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_subtract(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_subtract(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sub_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_subtract_inplace(A.size, &A[0], &B[0]) + clib.test_vm_double_subtract_inplace(int(A.size), &A[0], &B[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_mult(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_multiply(A.size, &A[0], 1, &B[0], &res[0]) + clib.test_vm_double_multiply(int(A.size), &A[0], 1, &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_mult_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_multiply_inplace(A.size, &A[0], 1, &B[0], 1) + clib.test_vm_double_multiply_inplace(int(A.size), &A[0], 1, &B[0], 1) @cython.boundscheck(False) @cython.wraparound(False) def vm_div(np.ndarray[double] A, np.ndarray[double] B): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_divide(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_divide(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_div_inplace(np.ndarray[double] A, np.ndarray[double] B): - clib.test_vm_double_divide_inplace(A.size, &A[0], &B[0]) + clib.test_vm_double_divide_inplace(int(A.size), &A[0], &B[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_cmult(np.ndarray[complex] A, np.ndarray[complex] B): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_double_complex_multiply(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_complex_multiply(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_cmult_conj(np.ndarray[complex] A, np.ndarray[complex] B): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_double_complex_conj_multiply(A.size, &A[0], &B[0], &res[0]) + clib.test_vm_double_complex_conj_multiply(int(A.size), &A[0], &B[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sq(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square(A.size, &A[0], &res[0]) + clib.test_vm_double_square(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sq_inplace(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square_inplace(A.size, &A[0]) + clib.test_vm_double_square_inplace(int(A.size), &A[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_sqrt(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square_root(A.size, &A[0], &res[0]) + clib.test_vm_double_square_root(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sqrt_inplace(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_square_root_inplace(A.size, &A[0]) + clib.test_vm_double_square_root_inplace(int(A.size), &A[0]) @cython.boundscheck(False) @cython.wraparound(False) def vm_cos(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_cosine(A.size, &A[0], &res[0]) + clib.test_vm_double_cosine(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_sin(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_sine(A.size, &A[0], &res[0]) + clib.test_vm_double_sine(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_cos_I_sin(np.ndarray[double] A): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_cosine_I_sine(A.size, &A[0], &res[0]) + clib.test_vm_cosine_I_sine(int(A.size), &A[0], &res[0]) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_exp(np.ndarray[double] A): cdef np.ndarray[double] res = np.zeros(A.size, dtype=float) - clib.test_vm_double_exp(A.size, &A[0], &res[0], 1) + clib.test_vm_double_exp(int(A.size), &A[0], &res[0], 1) return res @cython.boundscheck(False) @cython.wraparound(False) def vm_I_exp(np.ndarray[complex] A): cdef np.ndarray[complex] res = np.zeros(A.size, dtype=complex) - clib.test_vm_double_complex_exp(A.size, &A[0], &res[0]) + clib.test_vm_double_complex_exp(int(A.size), &A[0], &res[0]) return res diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index f349761d1..eccf9c7c9 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -395,13 +395,24 @@ def run( for index in method_index: method = self.methods[index] + + config_dict = self.config.get_int_dict() + # amp = core_simulator( + # method=method, + # spin_systems=self.spin_systems, + # transition_pathways=opt["precomputed_pathways"][index], + # transition_weights=opt["precomputed_weights"][index], + # **config_dict, + # **kwargs, + # ) + + # Chunk spin systems for multiple jobs spin_sys = get_chunks(self.spin_systems, n_jobs) - # Chunking transition pathways and weights form the optimization dictionary + # Chunk transition pathways and weights form the optimization dictionary pathways = get_chunks(opt["precomputed_pathways"][index], n_jobs) weights = get_chunks(opt["precomputed_weights"][index], n_jobs) - config_dict = self.config.get_int_dict() jobs = ( delayed(core_simulator)( method=method, @@ -418,6 +429,7 @@ def run( verbose=verbose, backend="loky", )(jobs) + amp = amp[0] B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density w_ref = method.channels[0].B0_to_ref_freq(B0) * 1e6 @@ -574,17 +586,8 @@ def _package_amp_after_simulation(self, method, amp, pack_as_csdm): bool pack_as_csdm: Packages the simulated spectrum as a CSDM object if true. Otherwise kept as a numpy array. """ - if isinstance(amp[0], list): - simulated_dataset = [] - for item in amp: - simulated_dataset += item - if isinstance(amp[0], np.ndarray): - simulated_dataset = [np.asarray(amp).sum(axis=0)] - method.simulation = ( - self._as_csdm_object(simulated_dataset, method) - if pack_as_csdm - else np.asarray(simulated_dataset) + self._as_csdm_object(amp, method) if pack_as_csdm else np.asarray(amp) ) diff --git a/tests/test_orientations_and_triangle_interpolation.py b/tests/test_orientations_and_triangle_interpolation.py index e19e5054c..97054de05 100644 --- a/tests/test_orientations_and_triangle_interpolation.py +++ b/tests/test_orientations_and_triangle_interpolation.py @@ -200,9 +200,8 @@ def test_triangle_interpolation(): # plt.title("1D interpolation, span(0, 100)") # plt.legend() # plt.tight_layout() - # plt.show() # plt.savefig(f"figs/fig_1D_{i}_{scl}.pdf") - # plt.figure().clear() + # plt.close() assert np.allclose(amp_py, amp_c.real, atol=1e-15)