diff --git a/LICENSE b/LICENSE index 77e93cc4e..018655a6e 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ BSD 3-Clause License -Copyright (c) 2019-2020, Mrsimulator Developors +Copyright (c) 2019-2021, Mrsimulator Developors All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/docs/configuring_simulator.rst b/docs/configuring_simulator.rst index f6618f7f3..a7047f853 100644 --- a/docs/configuring_simulator.rst +++ b/docs/configuring_simulator.rst @@ -51,9 +51,9 @@ of the config attributes is as follows, ... >>> sim = Simulator() >>> sim.config - ConfigSimulator(number_of_sidebands=64, integration_volume='octant', integration_density=70, decompose_spectrum='none') + ConfigSimulator(number_of_sidebands=64, integration_type='Alderman', integration_volume='octant', integration_density=70, decompose_spectrum='none') -Here, the configurable attributes are ``number_of_sidebands``, +Here, the configurable attributes are ``number_of_sidebands``, ``integration_type``, ``integration_volume``, ``integration_density``, and ``decompose_spectrum``. diff --git a/docs/credits/license.rst b/docs/credits/license.rst index 85f3a77d3..38a314eb3 100644 --- a/docs/credits/license.rst +++ b/docs/credits/license.rst @@ -6,7 +6,7 @@ Mrsimulator License Mrsimulator is licensed under BSD 3-Clause License -Copyright (c) 2019-2020, Mrsimulator Developors, +Copyright (c) 2019-2021, Mrsimulator Developors, All rights reserved. diff --git a/setup.py b/setup.py index 77517f533..4c82510f0 100644 --- a/setup.py +++ b/setup.py @@ -181,12 +181,13 @@ def message(lib): source = [ "src/c_lib/lib/angular_momentum.c", "src/c_lib/lib/interpolation.c", + "src/c_lib/lib/method.c", "src/c_lib/lib/mrsimulator.c", "src/c_lib/lib/octahedron.c", - "src/c_lib/lib/simulation.c", "src/c_lib/lib/powder_setup.c", "src/c_lib/lib/schemes.c", - "src/c_lib/lib/method.c", + "src/c_lib/lib/simulation.c", + "src/c_lib/lib/spherical.c", ] ext = ".pyx" if USE_CYTHON else ".c" diff --git a/src/c_lib/base/base_model.pxd b/src/c_lib/base/base_model.pxd index d37947dc5..18fb78117 100644 --- a/src/c_lib/base/base_model.pxd +++ b/src/c_lib/base/base_model.pxd @@ -2,7 +2,7 @@ # # nmr_method.pxd # -# @copyright Deepansh J. Srivastava, 2019-2020. +# @copyright Deepansh J. Srivastava, 2019-2021. # Created by Deepansh J. Srivastava. # Contact email = srivastava.89@osu.edu # @@ -17,9 +17,10 @@ cdef extern from "schemes.h": unsigned int total_orientations MRS_averaging_scheme *MRS_create_averaging_scheme( + unsigned int integration_type, unsigned int integration_density, - bool_t allow_fourth_rank, - unsigned int integration_volume) + unsigned int integration_volume, + bool_t allow_fourth_rank) void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) @@ -112,6 +113,7 @@ cdef extern from "simulation.h": # The transition as transition[0] = mi and transition[1] = mf float *transition, + unsigned int integration_type, # 0-octahedron, 1-sphere int integration_density, unsigned int integration_volume, # 0-octant, 1-hemisphere, 2-sphere bool_t interpolation, diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index 737f659f2..57f3ccd36 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -15,6 +15,7 @@ def one_d_spectrum(method, list spin_systems, int verbose=0, unsigned int number_of_sidebands=90, + unsigned int integration_type=1, unsigned int integration_density=72, unsigned int decompose_spectrum=0, unsigned int integration_volume=1, @@ -84,8 +85,8 @@ def one_d_spectrum(method, # create averaging scheme _____________________________________________________ cdef clib.MRS_averaging_scheme *the_averaging_scheme the_averaging_scheme = clib.MRS_create_averaging_scheme( - integration_density=integration_density, allow_fourth_rank=allow_fourth_rank, - integration_volume=integration_volume + integration_type=integration_type, integration_density=integration_density, + integration_volume=integration_volume, allow_fourth_rank=allow_fourth_rank ) # create spectral dimensions _______________________________________________ diff --git a/src/c_lib/include/angular_momentum.h b/src/c_lib/include/angular_momentum.h index 838f8d12a..79f149086 100644 --- a/src/c_lib/include/angular_momentum.h +++ b/src/c_lib/include/angular_momentum.h @@ -2,7 +2,7 @@ // // angular_momentum.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava. // Contact email = srivastava.89@osu.edu diff --git a/src/c_lib/include/array.h b/src/c_lib/include/array.h index 169b65fdc..dab6c2d93 100644 --- a/src/c_lib/include/array.h +++ b/src/c_lib/include/array.h @@ -2,7 +2,7 @@ // // array.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019 // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/config.h b/src/c_lib/include/config.h index 3ca2d72e2..7237d5f11 100644 --- a/src/c_lib/include/config.h +++ b/src/c_lib/include/config.h @@ -2,7 +2,7 @@ // // config.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Aug 10, 2019 // Contact email = srivastava.89@osu.edu // @@ -52,7 +52,8 @@ typedef float complex64[2]; #include "vm.h" // user definition -#define PI2 6.2831853072 +#define PI2 6.283185307179586 +#define PIby2 1.570796326794897 #define PI2I PI2 *I // #ifdef linux diff --git a/src/c_lib/include/frequency_component/spatial_orientation_tensor.h b/src/c_lib/include/frequency_component/spatial_orientation_tensor.h index 57a25654b..7ef37a53f 100644 --- a/src/c_lib/include/frequency_component/spatial_orientation_tensor.h +++ b/src/c_lib/include/frequency_component/spatial_orientation_tensor.h @@ -1,7 +1,7 @@ // // spatial_orientation_tensors.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Aug 26, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/frequency_component/spin_transition_function.h b/src/c_lib/include/frequency_component/spin_transition_function.h index 99be804ea..5d8a7c433 100644 --- a/src/c_lib/include/frequency_component/spin_transition_function.h +++ b/src/c_lib/include/frequency_component/spin_transition_function.h @@ -2,7 +2,7 @@ // // spin_transition_function.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/frequency_tensor.h b/src/c_lib/include/frequency_tensor.h index a96e7ce6e..f0546d0c5 100644 --- a/src/c_lib/include/frequency_tensor.h +++ b/src/c_lib/include/frequency_tensor.h @@ -2,7 +2,7 @@ // // frequency_component_function.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/interpolation.h b/src/c_lib/include/interpolation.h index 29762ca74..1264959e0 100644 --- a/src/c_lib/include/interpolation.h +++ b/src/c_lib/include/interpolation.h @@ -2,7 +2,7 @@ // // interpolation.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019 // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/isotopomer_ravel.h b/src/c_lib/include/isotopomer_ravel.h index 245ddf4de..6e2cc965b 100644 --- a/src/c_lib/include/isotopomer_ravel.h +++ b/src/c_lib/include/isotopomer_ravel.h @@ -2,7 +2,7 @@ // // isotopomer_ravel.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/method.h b/src/c_lib/include/method.h index 7179468ce..94d8dff73 100644 --- a/src/c_lib/include/method.h +++ b/src/c_lib/include/method.h @@ -3,7 +3,7 @@ // // method.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 15, 2020. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/mrsimulator.h b/src/c_lib/include/mrsimulator.h index 39f3c85cf..53f04930c 100644 --- a/src/c_lib/include/mrsimulator.h +++ b/src/c_lib/include/mrsimulator.h @@ -2,7 +2,7 @@ // // mrsimulator.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Jun 30, 2019 // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/octahedron.h b/src/c_lib/include/octahedron.h index ed5ac1c67..26a3792fd 100644 --- a/src/c_lib/include/octahedron.h +++ b/src/c_lib/include/octahedron.h @@ -2,7 +2,7 @@ // // powder_setup.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // @@ -10,15 +10,18 @@ #include "config.h" #include "interpolation.h" -extern void octahedronGetDirectionCosineSquareOverOctantAndWeights( - int nt, double *xr, double *yr, double *zr, double *amp); +extern void octahedronGetDirectionCosineSquareAndWeightsOverOctant( + int nt, double *restrict xr, double *restrict yr, double *restrict zr, + double *restrict amp); -extern void octahedronGetPolarAngleTrigOverAnOctant(int nt, double *cos_alpha, - double *cos_beta, - double *amp); +extern void octahedronGetPolarAngleTrigOverOctant(int nt, double *cos_alpha, + double *cos_beta, + double *amp); -extern void octahedronGetPolarAngleCosineAzimuthalAnglePhaseOverOctant( - int nt, void *exp_I_alpha, void *exp_I_beta, double *amp); +extern void octahedronGetComplexExpOfPolarAngleOverOctant(int nt, + void *exp_I_alpha, + void *exp_I_beta, + double *amp); extern void octahedronInterpolation(double *spec, double *freq, int nt, double *amp, int stride, int m); diff --git a/src/c_lib/include/powder_setup.h b/src/c_lib/include/powder_setup.h index 412c50ff2..87145d683 100644 --- a/src/c_lib/include/powder_setup.h +++ b/src/c_lib/include/powder_setup.h @@ -2,12 +2,12 @@ // // powder_setup.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // #include "config.h" -extern void octahedron_averaging_setup(int nt, void *exp_I_alpha, - void *exp_I_beta, double *amp); +extern void averaging_setup(unsigned int integration_type, int nt, + void *exp_I_alpha, void *exp_I_beta, double *amp); diff --git a/src/c_lib/include/schemes.h b/src/c_lib/include/schemes.h index 95de53280..1d9e50db3 100644 --- a/src/c_lib/include/schemes.h +++ b/src/c_lib/include/schemes.h @@ -2,7 +2,7 @@ // // schemes.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Sep 3, 2019. // Contact email = srivastava.89@osu.edu // @@ -40,8 +40,11 @@ typedef struct MRS_averaging_scheme { unsigned int total_orientations; /**< The total number of orientations. */ /** \privatesection */ - unsigned int integration_density; // number of triangles along the edge of - // the octahedron. + unsigned int integration_type; // Sampling type for integration + // 0-octaheron, 1-spherical + unsigned int integration_density; // number of triangles along the edge + // of the octahedron. + unsigned int integration_volume; // 0-octant, 1-hemisphere, 2-sphere. unsigned int octant_orientations; // number of unique orientations on the // face of an octant. double *amplitudes; // array of amplitude scaling per orientation. @@ -52,7 +55,6 @@ typedef struct MRS_averaging_scheme { double *wigner_4j_matrices; // wigner-d 4j matrix per orientation. // double *local_frequency; // buffer for local frequencies. // double *freq_offset; // buffer for local + sideband frequencies. - unsigned int integration_volume; // 0-octant, 1-hemisphere, 2-sphere. bool allow_fourth_rank; // If true, compute wigner matrices for wigner-d 4j. } MRS_averaging_scheme; @@ -77,8 +79,8 @@ typedef struct MRS_averaging_scheme { * @param integration_volume An enumeration. 0=octant, 1=hemisphere */ MRS_averaging_scheme *MRS_create_averaging_scheme( - unsigned int integration_density, bool allow_fourth_rank, - unsigned int integration_volume); + unsigned int integration_type, unsigned int integration_density, + unsigned int integration_volume, bool allow_fourth_rank); /** * Create a new orientation averaging scheme from given alpha and beta. @@ -105,6 +107,10 @@ MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( */ void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme); +static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, + complex128 *exp_I_beta, + bool allow_fourth_rank); + #endif // averaging_scheme_h #ifndef fftw_scheme_h diff --git a/src/c_lib/include/simulation.h b/src/c_lib/include/simulation.h index a0013dd62..633c6f10f 100644 --- a/src/c_lib/include/simulation.h +++ b/src/c_lib/include/simulation.h @@ -2,7 +2,7 @@ // // simulation.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // @@ -35,6 +35,7 @@ extern void mrsimulator_core( // powder orientation average // The number of triangle along the edge of octahedron. + unsigned int integration_type, // 0-octahedron, 1-sphere int integration_density, unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. bool interpolation, bool *freq_contrib, double *affine_matrix); diff --git a/src/c_lib/include/spherical.h b/src/c_lib/include/spherical.h new file mode 100644 index 000000000..17566bc52 --- /dev/null +++ b/src/c_lib/include/spherical.h @@ -0,0 +1,19 @@ +// -*- coding: utf-8 -*- +// +// spherical.h +// +// @copyright Deepansh J. Srivastava, 2019-2021. +// Created by Deepansh J. Srivastava, Feb 12, 2021. +// Contact email = srivastava.89@osu.edu +// + +#include "config.h" +#include "interpolation.h" + +extern void sphericalGetCoordinates(int nt, double *restrict x, + double *restrict y, double *restrict z); + +extern void sphericalGetComplexExpOfPolarAngleOverOctant(int nt, + void *exp_I_alpha, + void *exp_I_beta, + double *amp); diff --git a/src/c_lib/include/vm.h b/src/c_lib/include/vm.h index ef4b6fdde..c5bc95795 100644 --- a/src/c_lib/include/vm.h +++ b/src/c_lib/include/vm.h @@ -2,7 +2,7 @@ // // vm.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Jul 26, 2019. // Contact email = srivastava.89@osu.edu // @@ -243,15 +243,16 @@ static inline void vm_double_complex_multiply(int count, const void *restrict x, * res = cos(x) */ static inline void vm_double_cosine(int count, const double *restrict x, - double *restrict res) { + int stride_x, double *restrict res, + int stride_res) { // x = __builtin_assume_aligned(x, 32); // res = __builtin_assume_aligned(res, 32); while (count-- > 0) { - *res++ = cos(*x++); + *res = cos(*x); // res++; // x++; - // x += stride_x; - // res += stride_res; + x += stride_x; + res += stride_res; } } @@ -260,13 +261,14 @@ static inline void vm_double_cosine(int count, const double *restrict x, * res = sin(x) */ static inline void vm_double_sine(int count, const double *restrict x, - double *restrict res) { + int stride_x, double *restrict res, + int stride_res) { // x = __builtin_assume_aligned(x, 32); // res = __builtin_assume_aligned(res, 32); while (count-- > 0) { - *res++ = sin(*x++); - // x += stride_x; - // res += stride_res; + *res = sin(*x); + x += stride_x; + res += stride_res; } } diff --git a/src/c_lib/include/vm_common.h b/src/c_lib/include/vm_common.h index b84c65a8b..7a261fc64 100644 --- a/src/c_lib/include/vm_common.h +++ b/src/c_lib/include/vm_common.h @@ -2,7 +2,7 @@ // // vm_common.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Jul 26, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/include/vm_mkl.h b/src/c_lib/include/vm_mkl.h index b1c24476d..7d3f32f54 100644 --- a/src/c_lib/include/vm_mkl.h +++ b/src/c_lib/include/vm_mkl.h @@ -2,7 +2,7 @@ // // vm_mkl.h // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Jul 26, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/lib/interpolation.c b/src/c_lib/lib/interpolation.c index cd1ecfb09..0a6a55676 100644 --- a/src/c_lib/lib/interpolation.c +++ b/src/c_lib/lib/interpolation.c @@ -2,7 +2,7 @@ // // interpolation.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019. // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/lib/method.c b/src/c_lib/lib/method.c index 891f0d610..ec77cc915 100644 --- a/src/c_lib/lib/method.c +++ b/src/c_lib/lib/method.c @@ -2,7 +2,7 @@ // // method.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 15, 2020 // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/lib/mrsimulator.c b/src/c_lib/lib/mrsimulator.c index 7a400d1cb..0bc2d436a 100644 --- a/src/c_lib/lib/mrsimulator.c +++ b/src/c_lib/lib/mrsimulator.c @@ -2,7 +2,7 @@ // // mrsimulator.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Jun 9, 2019 // Contact email = srivastava.89@osu.edu // diff --git a/src/c_lib/lib/octahedron.c b/src/c_lib/lib/octahedron.c index dba370688..d7a424c79 100644 --- a/src/c_lib/lib/octahedron.c +++ b/src/c_lib/lib/octahedron.c @@ -2,18 +2,17 @@ // // octahedron.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019 // Contact email = srivastava.89@osu.edu // #include "octahedron.h" -void octahedronGetDirectionCosineSquareOverOctantAndWeights(int nt, double *xr, - double *yr, - double *zr, - double *amp) { - int i, j, k = 0; +void octahedronGetDirectionCosineSquareAndWeightsOverOctant( + int nt, double *restrict xr, double *restrict yr, double *restrict zr, + double *restrict amp) { + int i, j; double x2, y2, z2, r2, scale = (double)nt; /* Do the (x + y + z = nt) face of the octahedron @@ -28,23 +27,22 @@ void octahedronGetDirectionCosineSquareOverOctantAndWeights(int nt, double *xr, y2 = pow(i, 2); z2 = pow(j, 2); r2 = x2 + y2 + z2; - xr[k] = x2 / r2; - yr[k] = y2 / r2; - zr[k] = z2 / r2; - amp[k] = scale / (r2 * sqrt(r2)); - k++; + *xr++ = x2 / r2; + *yr++ = y2 / r2; + *zr++ = z2 / r2; + *amp++ = scale / (r2 * sqrt(r2)); } } - xr[k] = 0.0; - yr[k] = 0.0; - zr[k] = 1.0; + *xr = 0.0; + *yr = 0.0; + *zr = 1.0; r2 = (double)nt; - amp[k] = scale / (r2 * r2 * r2); + *amp = 1.0 / (r2 * r2); } -void octahedronGetPolarAngleTrigOverAnOctant(int nt, double *cos_alpha, - double *cos_beta, double *amp) { +void octahedronGetPolarAngleTrigOverOctant(int nt, double *cos_alpha, + double *cos_beta, double *amp) { int points = (nt + 1) * (nt + 2) / 2; double *xr = malloc_double(points); double *yr = malloc_double(points); @@ -53,7 +51,7 @@ void octahedronGetPolarAngleTrigOverAnOctant(int nt, double *cos_alpha, // The values xr = x^2, yr = y^2, zr = z^2, where x, y, and z are the // direction cosines. - octahedronGetDirectionCosineSquareOverOctantAndWeights(nt, xr, yr, zr, amp); + octahedronGetDirectionCosineSquareAndWeightsOverOctant(nt, xr, yr, zr, amp); // Evaluate sqrt of zr to get cos(beta) ---> cos(beta) = sqrt(z^2) vm_double_square_root(points, zr, cos_beta); @@ -83,14 +81,15 @@ void octahedronGetPolarAngleTrigOverAnOctant(int nt, double *cos_alpha, free(sin_beta); } -void octahedronGetPolarAngleCosineAzimuthalAnglePhaseOverOctant( - int nt, void *exp_I_alpha, void *exp_I_beta, double *amp) { +void octahedronGetComplexExpOfPolarAngleOverOctant(int nt, void *exp_I_alpha, + void *exp_I_beta, + double *amp) { int points = (nt + 1) * (nt + 2) / 2; double *xr = malloc_double(points); double *yr = malloc_double(points); double *zr = malloc_double(points); - octahedronGetDirectionCosineSquareOverOctantAndWeights(nt, xr, yr, zr, amp); + octahedronGetDirectionCosineSquareAndWeightsOverOctant(nt, xr, yr, zr, amp); // At this point the variables // xr = x^2, // yr = y^2, @@ -99,29 +98,29 @@ void octahedronGetPolarAngleCosineAzimuthalAnglePhaseOverOctant( // Cos beta .............................................................. // // cos(beta) = sqrt(z^2) ... (1) - // In terms of the variables, Eq (1) is given as sqrt(zr) - // Evaluate zr = sqrt(zr) ==>> sqrt(z^2) + // In terms of the variables, Eq (1) is given as sqrt(zr). + // Evaluate zr = sqrt(zr) ==>> sqrt(z^2). vm_double_square_root(points, zr, zr); - // Copy cos(beta), aka zr, to the even addresses of exp_I_beta + // Copy cos(beta), aka zr, to the even addresses of exp_I_beta. cblas_dcopy(points, zr, 1, (double *)exp_I_beta, 2); // Sin beta .............................................................. // // sin^2(beta) = x^2 + y^2 ... (2) - // In terms of the variables, Eq (2) is given as xr + yr + // In terms of the variables, Eq (2) is given as xr + yr. // Evaluate zr = xr + yr. vm_double_add(points, xr, yr, zr); // Take the square root of zr to get sin(beta). - // Evaluate zr = sqrt(zr) ==>> sqrt(x^2 + y^2) + // Evaluate zr = sqrt(zr) ==>> sqrt(x^2 + y^2). vm_double_square_root(points, zr, zr); - // Copy sin(beta), aka zr, to the odd addresses of exp_I_beta + // Copy sin(beta), aka zr, to the odd addresses of exp_I_beta. cblas_dcopy(points, zr, 1, (double *)exp_I_beta + 1, 2); // Evaluate squate root of xr - // xr = sqrt(xr) ==>> sqrt(x^2) + // xr = sqrt(xr) ==>> sqrt(x^2). vm_double_square_root(points, xr, xr); // Evaluate squate root of yr - // yr = sqrt(yr) ==>> sqrt(y^2) + // yr = sqrt(yr) ==>> sqrt(y^2). vm_double_square_root(points, yr, yr); // .. note diff --git a/src/c_lib/lib/powder_setup.c b/src/c_lib/lib/powder_setup.c index f8cfd5714..5d8cdfed0 100644 --- a/src/c_lib/lib/powder_setup.c +++ b/src/c_lib/lib/powder_setup.c @@ -2,7 +2,7 @@ // // powder_setup.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019 // Contact email = srivastava.89@osu.edu // @@ -10,10 +10,18 @@ #include "powder_setup.h" #include "octahedron.h" +#include "spherical.h" -void octahedron_averaging_setup(int nt, void *exp_I_alpha, void *exp_I_beta, - double *amp) { - // octahedronGetPolarAngleTrigOverAnOctant(nt, cos_alpha, cos_beta, amp); - octahedronGetPolarAngleCosineAzimuthalAnglePhaseOverOctant(nt, exp_I_alpha, - exp_I_beta, amp); +void averaging_setup(unsigned int integration_type, int nt, void *exp_I_alpha, + void *exp_I_beta, double *amp) { + // octahedronGetPolarAngleTrigOverOctant(nt, cos_alpha, cos_beta, amp); + if (integration_type == 0) { + octahedronGetComplexExpOfPolarAngleOverOctant(nt, exp_I_alpha, exp_I_beta, + amp); + return; + } + if (integration_type == 1) { + sphericalGetComplexExpOfPolarAngleOverOctant(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 1c821fc57..796872a65 100644 --- a/src/c_lib/lib/schemes.c +++ b/src/c_lib/lib/schemes.c @@ -2,13 +2,77 @@ // // schemes.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava. // Contact email = srivastava.89@osu.edu // #include "schemes.h" +/* Free the memory from the mrsimulator plan associated with the spherical + * averaging scheme */ +void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) { + free(scheme->amplitudes); + free(scheme->exp_Im_alpha); + free(scheme->w2); + free(scheme->w4); + free(scheme->wigner_2j_matrices); + free(scheme->wigner_4j_matrices); +} + +/* Create a new orientation averaging scheme. */ +MRS_averaging_scheme *MRS_create_averaging_scheme( + unsigned int integration_type, unsigned int integration_density, + unsigned int integration_volume, bool allow_fourth_rank) { + MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); + + scheme->integration_density = integration_density; + scheme->integration_volume = integration_volume; + scheme->allow_fourth_rank = allow_fourth_rank; + + scheme->octant_orientations = + ((integration_density + 1) * (integration_density + 2)) / 2; + + /* 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); + complex128 *exp_I_beta = malloc_complex128(scheme->octant_orientations); + scheme->amplitudes = malloc_double(scheme->octant_orientations); + + averaging_setup(integration_type, integration_density, + &scheme->exp_Im_alpha[3 * scheme->octant_orientations], + exp_I_beta, scheme->amplitudes); + + averaging_scheme_setup(scheme, exp_I_beta, allow_fourth_rank); + return scheme; +} + +/* 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_fourth_rank) { + MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); + + scheme->octant_orientations = n_angles; + scheme->integration_volume = 0; + scheme->total_orientations = n_angles; + + scheme->exp_Im_alpha = malloc_complex128(4 * scheme->total_orientations); + complex128 *exp_I_beta = malloc_complex128(scheme->total_orientations); + scheme->amplitudes = weight; + + /* Calculate cos(α) + isin(α) from α. .................................... */ + vm_cosine_I_sine(n_angles, alpha, + scheme->exp_Im_alpha[3 * scheme->octant_orientations]); + + /* Calculate cos(β) + isin(β) from β. .................................... */ + vm_cosine_I_sine(n_angles, beta, exp_I_beta); + + averaging_scheme_setup(scheme, exp_I_beta, allow_fourth_rank); + return scheme; +} + static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, complex128 *exp_I_beta, bool allow_fourth_rank) { @@ -109,71 +173,6 @@ static inline void averaging_scheme_setup(MRS_averaging_scheme *scheme, } } -/* Free the memory from the mrsimulator plan associated with the spherical - * averaging scheme */ -void MRS_free_averaging_scheme(MRS_averaging_scheme *scheme) { - free(scheme->amplitudes); - free(scheme->exp_Im_alpha); - free(scheme->w2); - free(scheme->w4); - free(scheme->wigner_2j_matrices); - free(scheme->wigner_4j_matrices); -} - -/* Create a new orientation averaging scheme. */ -MRS_averaging_scheme *MRS_create_averaging_scheme( - unsigned int integration_density, bool allow_fourth_rank, - unsigned int integration_volume) { - MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); - - scheme->integration_density = integration_density; - scheme->integration_volume = integration_volume; - scheme->allow_fourth_rank = allow_fourth_rank; - - scheme->octant_orientations = - ((integration_density + 1) * (integration_density + 2)) / 2; - - /* 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); - complex128 *exp_I_beta = malloc_complex128(scheme->octant_orientations); - scheme->amplitudes = malloc_double(scheme->octant_orientations); - - octahedron_averaging_setup( - integration_density, - &scheme->exp_Im_alpha[3 * scheme->octant_orientations], exp_I_beta, - scheme->amplitudes); - - averaging_scheme_setup(scheme, exp_I_beta, allow_fourth_rank); - return scheme; -} - -/* 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_fourth_rank) { - MRS_averaging_scheme *scheme = malloc(sizeof(MRS_averaging_scheme)); - - scheme->octant_orientations = n_angles; - scheme->integration_volume = 0; - scheme->total_orientations = n_angles; - - scheme->exp_Im_alpha = malloc_complex128(4 * scheme->total_orientations); - complex128 *exp_I_beta = malloc_complex128(scheme->total_orientations); - scheme->amplitudes = weight; - - /* Calculate cos(α) + isin(α) from α. .................................... */ - vm_cosine_I_sine(n_angles, alpha, - scheme->exp_Im_alpha[3 * scheme->octant_orientations]); - - /* Calculate cos(β) + isin(β) from β. .................................... */ - vm_cosine_I_sine(n_angles, beta, exp_I_beta); - - averaging_scheme_setup(scheme, exp_I_beta, allow_fourth_rank); - return scheme; -} - /* ----------------------------------------------------------------------- */ /* fftw routine setup .................................................... */ /* ....................................................................... */ diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index 3947fc74c..4dac9d135 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -2,7 +2,7 @@ // // sideband_simulator.c // -// @copyright Deepansh J. Srivastava, 2019-2020. +// @copyright Deepansh J. Srivastava, 2019-2021. // Created by Deepansh J. Srivastava, Apr 11, 2019 // Contact email = srivastava.89@osu.edu // @@ -424,6 +424,8 @@ void mrsimulator_core( float *transition, // powder orientation average + unsigned int integration_type, // Sampling type for integration + // 0-octaheron, 1-spherical int integration_density, // The number of triangle along the edge // of octahedron unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. @@ -448,8 +450,9 @@ void mrsimulator_core( number_of_sidebands = 1; } - MRS_averaging_scheme *scheme = MRS_create_averaging_scheme( - integration_density, allow_fourth_rank, integration_volume); + MRS_averaging_scheme *scheme = + MRS_create_averaging_scheme(integration_type, integration_density, + integration_volume, allow_fourth_rank); MRS_fftw_scheme *fftw_scheme = create_fftw_scheme(scheme->total_orientations, number_of_sidebands); diff --git a/src/c_lib/lib/spherical.c b/src/c_lib/lib/spherical.c new file mode 100644 index 000000000..c253019e8 --- /dev/null +++ b/src/c_lib/lib/spherical.c @@ -0,0 +1,101 @@ +// -*- coding: utf-8 -*- +// +// spherical.c +// +// @copyright Deepansh J. Srivastava, 2019-2021. +// Created by Deepansh J. Srivastava, Feb 12, 2021 +// Contact email = srivastava.89@osu.edu +// + +#include "spherical.h" + +void sphericalGetCoordinates(int nt, double *restrict x, double *restrict y, + double *restrict z) { + int i, j; + double temp1, temp2, scale = (double)nt, x_o, y_o, z_o; + + for (j = 0; j <= nt - 1; j++) { + for (i = 0; i <= nt - j; i++) { + z_o = (double)j / scale; + y_o = (double)i / scale; + x_o = (double)(nt - i - j) / scale; + temp1 = PIby2 * (y_o / (y_o + x_o)); + *z = 2.0 * z_o * (1.0 - 0.5 * z_o); + temp2 = sqrt(1 - *z * *z); + *x++ = temp2 * cos(temp1); + *y++ = temp2 * sin(temp1); + z++; + } + } + *x = 0.0; + *y = 0.0; + *z = 1.0; +} + +void sphericalGetComplexExpOfPolarAngleOverOctant(int nt, void *exp_I_alpha, + void *exp_I_beta, + double *amp) { + /** Uniform area trioangle coordinates over the Octahedron. + * + * z_o = (double)j / scale (1a) + * y_o = (double)i / scale (1b) + * x_o = (double)(nt-i-j) / scale (1c) + * + * Area preservong mapping of the (x_o, y_o, z_o) coordinates over the + * octahedron to (x, y, z) over the sphere. + * + * z = 2.0 * z_o * (1.0 - 0.5 * z_o) (2a) + * x = sqrt(1 - z^2) cos( (π/2) (y_o / (y_o + x_o)) ) (2b) + * y = sqrt(1 - z^2) sin( (π/2) (y_o / (y_o + x_o)) ) (2c) + * + * Cosines of alpha and beta from the (x, y, z) coordinates + * + * cos_beta = z (3a) + * sin_beta = sqrt(1 - z^2) (3b) + * cos_alpha = sqrt(x^2 / (x^2 + y^2)) (3c) + * sin_alpha = sqrt(y^2 / (x^2 + y^2)) (3d) + * + * From Eqs. (2b), (2c), and (3b) + * x^2 + y^2 = 1 - z^2 = (sin_beta)^2 (4) + * + * Substituting Eqs. (4) and (2b) into (3c) gives + * cos_alpha = x / sin_beta (5a) + * cos_alpha = sin_beta cos( arg ) / sin_beta (5b) + * cos_alpha = cos( arg ) (5b) + * + * where, + * arg = (π/2) (y_o / (y_o + x_o)) (6a) + * = (π/2) ( i / (nt - j) ) (6b) + * + * Similarly, + * sin_alpha = y / sin_beta (7a) + * sin_alpha = sin_beta sin( arg ) / sin_beta (7b) + * sin_alpha = sin( arg ) (7c) + **/ + + int i, j, points = (nt + 1) * (nt + 2) / 2, index; + double temp, arg, z, scale = (double)nt; + double *exp_I_alpha_ = (double *)exp_I_alpha; + double *exp_I_beta_ = (double *)exp_I_beta; + + vm_double_ones(points, amp); + + for (j = 0; j <= nt - 1; j++) { + for (i = 0; i <= nt - j; i++) { + z = (double)j / scale; // z-coordinate over octahedron + temp = (double)(nt - j); // x + y + + z = 2.0 * z * (1.0 - 0.5 * z); // z-coordinate over sphere + *exp_I_beta_++ = z; + *exp_I_beta_++ = sqrt(1.0 - z * z); + + arg = PIby2 * (double)i / temp; + *exp_I_alpha_++ = cos(arg); + *exp_I_alpha_++ = sin(arg); + } + } + *exp_I_beta_++ = 1.0; + *exp_I_beta_ = 0.0; + *exp_I_alpha_++ = 1.0; + *exp_I_alpha_ = 0.0; +} diff --git a/src/c_lib/sandbox/sandbox.pxd b/src/c_lib/sandbox/sandbox.pxd index 767a2139e..145876516 100644 --- a/src/c_lib/sandbox/sandbox.pxd +++ b/src/c_lib/sandbox/sandbox.pxd @@ -2,7 +2,7 @@ # # sandbox.pxd # -# @copyright Deepansh J. Srivastava, 2019-2020. +# @copyright Deepansh J. Srivastava, 2019-2021. # Created by Deepansh J. Srivastava # Contact email = srivastava.89@osu.edu # @@ -12,13 +12,15 @@ from libcpp cimport bool as bool_t cdef extern from "schemes.h": ctypedef struct MRS_averaging_scheme: unsigned int total_orientations + unsigned int integration_type unsigned int integration_density unsigned int integration_volume MRS_averaging_scheme * MRS_create_averaging_scheme( + unsigned int integration_type, unsigned int integration_density, - bool_t allow_fourth_rank, - unsigned int integration_volume) + unsigned int integration_volume, + bool_t allow_fourth_rank) MRS_averaging_scheme *MRS_create_averaging_scheme_from_alpha_beta( double *alpha, double *beta, diff --git a/src/c_lib/sandbox/sandbox.pyx b/src/c_lib/sandbox/sandbox.pyx index 262cd07ef..466c10302 100644 --- a/src/c_lib/sandbox/sandbox.pyx +++ b/src/c_lib/sandbox/sandbox.pyx @@ -69,7 +69,7 @@ cdef class AveragingScheme: cdef clib.MRS_averaging_scheme *scheme cdef bool_t allow_fourth_rank - def __init__(self, int integration_density, integration_volume='octant', bool_t allow_fourth_rank=False): + def __init__(self, int integration_density, integration_type='sphere', integration_volume='octant', bool_t allow_fourth_rank=False): """Create the octahedral interpolation averaging scheme for the simulation. Args: @@ -79,9 +79,9 @@ cdef class AveragingScheme: """ self.allow_fourth_rank = allow_fourth_rank integration_volume_ = 0 - if integration_volume == 'hemisphere': - integration_volume_=1 - self.scheme = clib.MRS_create_averaging_scheme(integration_density, + integration_volume_ = 1 if integration_volume == 'hemisphere' else 0 + integration_type_ = 1 if integration_type == 'sphere' else 0 + self.scheme = clib.MRS_create_averaging_scheme(integration_type, integration_density, allow_fourth_rank, integration_volume_) @property @@ -101,9 +101,9 @@ cdef class AveragingScheme: def integration_density(self, value): if isinstance(value, int): if value > 0: - self.scheme = clib.MRS_create_averaging_scheme(value, - self.allow_fourth_rank, - self.scheme.integration_volume) + self.scheme = clib.MRS_create_averaging_scheme(0, value, + self.scheme.integration_volume, + self.allow_fourth_rank) return raise ValueError(f"Expecting a positive integer, found {value}.") @@ -116,9 +116,9 @@ cdef class AveragingScheme: def integration_volume(self, value): if value in __integration_volume_enum__.keys(): self.scheme.integration_volume = __integration_volume_enum__[value] - self.scheme = clib.MRS_create_averaging_scheme( - self.scheme.integration_density, - self.allow_fourth_rank, value + self.scheme = clib.MRS_create_averaging_scheme(0, + self.scheme.integration_density, value, + self.allow_fourth_rank ) return raise ValueError( diff --git a/src/c_lib/test/test.pxd b/src/c_lib/test/test.pxd index 7e40843ca..51da9d8d9 100644 --- a/src/c_lib/test/test.pxd +++ b/src/c_lib/test/test.pxd @@ -2,7 +2,7 @@ # # test.pxd # -# @copyright Deepansh J. Srivastava, 2019-2020. +# @copyright Deepansh J. Srivastava, 2019-2021. # Created by Deepansh J. Srivastava. # Contact email = srivastava.89@osu.edu # @@ -41,7 +41,8 @@ cdef extern from "angular_momentum.h": cdef extern from "powder_setup.h": - void octahedron_averaging_setup( + void averaging_setup( + unsigned int integration_type, int nt, double complex *exp_I_alpha, double complex *exp_I_beta, @@ -157,6 +158,7 @@ cdef extern from "simulation.h": # The transition as transition[0] = mi and transition[1] = mf double *transition, + unsigned int integration_type, # 0-octahedron, 1-sphere int integration_density, unsigned int integration_volume, # 0-octant, 1-hemisphere, 2-sphere. bool_t interpolation diff --git a/src/c_lib/test/test.pyx b/src/c_lib/test/test.pyx index 0a43e66e2..4ce589d93 100644 --- a/src/c_lib/test/test.pyx +++ b/src/c_lib/test/test.pyx @@ -157,7 +157,8 @@ def pre_phase_components(unsigned int number_of_sidebands, double sample_rotatio @cython.boundscheck(False) @cython.wraparound(False) -def cosine_of_polar_angles_and_amplitudes(int integration_density=72): +def cosine_of_polar_angles_and_amplitudes(int integration_density=72, + unsigned int integration_type=1): r""" Calculate the direction cosines and the related amplitudes for the positive quadrant of the sphere. The direction cosines corresponds to @@ -191,7 +192,9 @@ def cosine_of_polar_angles_and_amplitudes(int integration_density=72): cdef np.ndarray[double complex] exp_I_beta = np.empty(octant_orientations, dtype=np.complex128) cdef np.ndarray[double] amp = np.empty(octant_orientations, dtype=np.float64) - clib.octahedron_averaging_setup(nt, &exp_I_alpha[0], &exp_I_beta[0], &[0]) + clib.averaging_setup( + integration_type, nt, &exp_I_alpha[0], &exp_I_beta[0], &[0] + ) return exp_I_alpha, exp_I_beta, amp diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index 6e857549e..7594d1cad 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -288,6 +288,7 @@ def json(self, include_methods: bool = False, include_version: bool = False): >>> pprint(sim.json()) {'config': {'decompose_spectrum': 'none', 'integration_density': 70, + 'integration_type': 'Alderman', 'integration_volume': 'octant', 'number_of_sidebands': 64}, 'spin_systems': [{'abundance': '100 %', diff --git a/src/mrsimulator/simulator/config.py b/src/mrsimulator/simulator/config.py index 409d2bdb8..f22f84de1 100644 --- a/src/mrsimulator/simulator/config.py +++ b/src/mrsimulator/simulator/config.py @@ -15,6 +15,8 @@ __integration_volume_enum__ = {"octant": 0, "hemisphere": 1} __integration_volume_octants__ = [1, 4] +__integration_type_enum__ = {"Alderman": 0, "spherical": 1} + class ConfigSimulator(BaseModel): r""" @@ -28,23 +30,27 @@ class ConfigSimulator(BaseModel): simulation. The value cannot be zero or negative. The default value is 64. + integration_type: enum (optional). + Orientation sampling type. The valid literals of this enumeration are + + - ``Spherical`` (default), and + - ``Aldermans`` + integration_volume: enum (optional). - The value is the volume over which the solid-state spectral frequency - integration is performed. The valid literals of this enumeration are + The volume over which the solid-state spectral frequency integral is evaluated. + The valid literals of this enumeration are - ``octant`` (default), and - ``hemisphere`` integration_density: int (optional). - The value represents the integration density or equivalently the number of - orientations over which the frequency integration is performed within a given - volume. If :math:`n` is the integration_density, then the total number of - orientation is given as + The integration density or equivalently the number of orientations sampled over + an octant. If :math:`n` is the integration_density, then the total number of + orientation, :math:`N`, sampled over an octant is .. math:: - n_\text{octants} \frac{(n+1)(n+2)}{2}, + N = \frac{(n+1)(n+2)}{2}. - where :math:`n_\text{octants}` is the number of octants in the given volume. The default value is 70. decompose_spectrum: enum (optional). @@ -62,12 +68,14 @@ class ConfigSimulator(BaseModel): >>> a = Simulator() >>> a.config.number_of_sidebands = 128 + >>> a.config.integration_type = 'Alderman' >>> a.config.integration_density = 96 >>> a.config.integration_volume = 'hemisphere' >>> a.config.decompose_spectrum = 'spin_system' """ number_of_sidebands: int = Field(default=64, gt=0) + integration_type: Literal["Alderman", "spherical"] = "Alderman" integration_volume: Literal["octant", "hemisphere"] = "octant" integration_density: int = Field(default=70, gt=0) decompose_spectrum: Literal["none", "spin_system"] = "none" @@ -83,6 +91,7 @@ def get_int_dict(self): py_dict["decompose_spectrum"] = __decompose_spectrum_enum__[ self.decompose_spectrum ] + py_dict["integration_type"] = __integration_type_enum__[self.integration_type] return py_dict # averaging scheme. This contains the c pointer used in frequency evaluation diff --git a/src/mrsimulator/simulator/tests/test_config.py b/src/mrsimulator/simulator/tests/test_config.py index 75645ca3e..7fc0f0289 100644 --- a/src/mrsimulator/simulator/tests/test_config.py +++ b/src/mrsimulator/simulator/tests/test_config.py @@ -12,6 +12,21 @@ def test_config(): # set config b = Simulator() + assert b.config.dict() == { + "decompose_spectrum": "none", + "number_of_sidebands": 64, + "integration_type": "Alderman", + "integration_volume": "octant", + "integration_density": 70, + } + assert b.config.get_int_dict() == { + "decompose_spectrum": 0, + "number_of_sidebands": 64, + "integration_type": 0, + "integration_volume": 0, + "integration_density": 70, + } + error = "value is not a valid dict" with pytest.raises(ValueError, match=f".*{error}.*"): b.config = 5 @@ -33,6 +48,15 @@ def test_config(): with pytest.raises(ValueError, match=f".*{error}.*"): a.config.number_of_sidebands = 0 + # integration type + assert a.config.integration_type == "Alderman" + a.config.integration_type = "spherical" + assert a.config.integration_type == "spherical" + + error = "unexpected value; permitted: 'Alderman', 'spherical'" + with pytest.raises(ValueError, match=f".*{error}.*"): + a.config.integration_type = "haha" + # integration density assert a.config.integration_density == 70 a.config.integration_density = 20 @@ -68,6 +92,7 @@ def test_config(): assert a.config.dict() == { "decompose_spectrum": "spin_system", "number_of_sidebands": 10, + "integration_type": "spherical", "integration_volume": "hemisphere", "integration_density": 20, } @@ -75,6 +100,7 @@ def test_config(): assert a.config.get_int_dict() == { "decompose_spectrum": 1, "number_of_sidebands": 10, + "integration_type": 1, "integration_volume": 1, "integration_density": 20, } diff --git a/src/mrsimulator/simulator/tests/test_simulator.py b/src/mrsimulator/simulator/tests/test_simulator.py index f7b9083b6..4b144d6f7 100644 --- a/src/mrsimulator/simulator/tests/test_simulator.py +++ b/src/mrsimulator/simulator/tests/test_simulator.py @@ -46,6 +46,7 @@ def test_equality(): "spin_systems": [{"abundance": "100 %", "sites": []}], "config": { "decompose_spectrum": "none", + "integration_type": "Alderman", "integration_density": 70, "integration_volume": "octant", "number_of_sidebands": 64, @@ -59,6 +60,7 @@ def test_equality(): "methods": [], "config": { "number_of_sidebands": 64, + "integration_type": "Alderman", "integration_volume": "octant", "integration_density": 70, "decompose_spectrum": "none", @@ -138,6 +140,7 @@ def test_simulator_1(): ], "config": { "decompose_spectrum": "none", + "integration_type": "Alderman", "integration_density": 70, "integration_volume": "octant", "number_of_sidebands": 64, diff --git a/tests/python_test_for_c_code/orientation.py b/tests/python_test_for_c_code/orientation.py index d2e726888..f0c00f62b 100644 --- a/tests/python_test_for_c_code/orientation.py +++ b/tests/python_test_for_c_code/orientation.py @@ -5,101 +5,161 @@ __email__ = "srivastava.89@osu.edu" -# def octahedral_coordinate(nt: int): - -# # Do the (x + y + z = nt) face of the octahedron -# # z -> 0 to nt-1 -# # y -> 0 to nt-z -# # x -> nt - y - z - -# n = int((nt + 1) * (nt + 2) / 2) -# x = np.empty(n, dtype=np.float64) -# y = np.empty(n, dtype=np.float64) -# z = np.empty(n, dtype=np.float64) - -# k = 0 -# for j in range(nt + 1): -# for i in range(nt - j + 1): -# # x = nt-i-j; -# # y = i; -# # z = j; -# x[k] = nt - i - j -# y[k] = i -# z[k] = j -# k += 1 - -# return x, y, z - - -def octahedral_direction_cosine_squares_and_amplitudes(nt: int): - - # Do the (x + y + z = nt) face of the octahedron - # z -> 0 to nt-1 - # y -> 0 to nt-z - # x -> nt - y - z - - n = int((nt + 1) * (nt + 2) / 2) - x = np.empty(n, dtype=np.float64) - y = np.empty(n, dtype=np.float64) - z = np.empty(n, dtype=np.float64) - - k = 0 - for j in range(nt + 1): - for i in range(nt - j + 1): - # x = nt-i-j; - # y = i; - # z = j; - x[k] = nt - i - j - y[k] = i - z[k] = j - k += 1 - - x *= x - y *= y - z *= z - - r2 = x + y + z - - x /= r2 - y /= r2 - z /= r2 - amp = nt / (r2 * np.sqrt(r2)) - - return x, y, z, amp - - -def cosine_of_polar_angles_and_amplitudes(integration_density: int = 72): - r""" - Calculates and return the direction cosines and the related amplitudes for - the positive quadrant of the sphere. The direction cosines corresponds to - angle $\alpha$ and $\beta$, where $\alpha$ is the azimuthal angle and - $\beta$ is the polar angle. The amplitudes are evaluated as $\frac{1}{r^3}$ - where $r$ is the distance from the origin to the face of the unit - octahedron in the positive quadrant along the line given by the values of - $\alpha$ and $\beta$. - - :ivar integration_density: - The value is an integer which represents the frequency of class I - geodesic polyhedra. These polyhedra are used in calculating the - spherical average. Presently we only use octahedral as the frequency1 - polyhedra. As the frequency of the geodesic polyhedron increases, the - polyhedra approach a sphere geometry. A higher frequency will result in - a better powder averaging. The default value is 72. - Read more on the `Geodesic polyhedron - `_. - - :return cos_alpha: The cosine of the azimuthal angle. - :return cos_beta: The cosine of the polar angle. - :return amp: The amplitude at the given $\alpha$ and $\beta$. - """ - nt = integration_density - xr, yr, zr, amp = octahedral_direction_cosine_squares_and_amplitudes(nt) - - cos_beta = np.sqrt(zr) - cos_alpha = np.zeros(xr.size, dtype=np.float64) - cos_alpha[:-1] = np.sqrt(xr[:-1] / (xr[:-1] + yr[:-1])) - cos_alpha[-1] = 0.0 - return cos_alpha, cos_beta, amp +class Octahedron: + def __init__(self, nt): + """Octahedron class. nt is the number of triangles along the edge of the + octahedron.""" + self.nt = nt + + def coordinates(self, nt: int): + """Coordinates over the position octant face of a octahedron. + + Args: + nt: Number of triangles along the edge of the octahedron face. + """ + # Do the (x + y + z = nt) face of the octahedron + # z -> 0 to nt-1 + # y -> 0 to nt-z + # x -> nt - y - z + + n = int((nt + 1) * (nt + 2) / 2) + x = np.empty(n, dtype=np.float64) + y = np.empty(n, dtype=np.float64) + z = np.empty(n, dtype=np.float64) + + k = 0 + for j in range(nt): + for i in range(nt - j + 1): + x[k] = nt - i - j + y[k] = i + z[k] = j + k += 1 + + x[k] = 0 + y[k] = 0 + z[k] = nt + + return x, y, z + + def projected_coordinates(self, nt: int): + """The coordinates over the face of octahedron projected onto the sphere.""" + x, y, z = self.coordinates(nt) + r = np.sqrt(x * x + y * y + z * z) + return x / r, y / r, z / r + + def _direction_cosine_squares_and_amplitudes(self, nt: int): + x, y, z = self.coordinates(nt) + + x *= x + y *= y + z *= z + + r2 = x + y + z + + x /= r2 + y /= r2 + z /= r2 + amp = nt / (r2 * np.sqrt(r2)) + + return x, y, z, amp + + def cosines_and_amplitudes(self, integration_density: int = 72): + r""" + Calculates and return the direction cosines and the related amplitudes for + the positive quadrant of the sphere. The direction cosines corresponds to + angle $\alpha$ and $\beta$, where $\alpha$ is the azimuthal angle and + $\beta$ is the polar angle. The amplitudes are evaluated as $\frac{1}{r^3}$ + where $r$ is the distance from the origin to the face of the unit + octahedron in the positive quadrant along the line given by the values of + $\alpha$ and $\beta$. + + :ivar integration_density: + The value is an integer which represents the frequency of class I + geodesic polyhedra. These polyhedra are used in calculating the + spherical average. Presently we only use octahedral as the frequency1 + polyhedra. As the frequency of the geodesic polyhedron increases, the + polyhedra approach a sphere geometry. A higher frequency will result in + a better powder averaging. The default value is 72. + Read more on the `Geodesic polyhedron + `_. + + :return cos_alpha: The cosine of the azimuthal angle. + :return cos_beta: The cosine of the polar angle. + :return amp: The amplitude at the given $\alpha$ and $\beta$. + """ + nt = integration_density + xr, yr, zr, amp = self._direction_cosine_squares_and_amplitudes(nt) + + cos_beta = np.sqrt(zr) + cos_alpha = np.zeros(xr.size, dtype=np.float64) + cos_alpha[:-1] = np.sqrt(xr[:-1] / (xr[:-1] + yr[:-1])) + cos_alpha[-1] = 0.0 + return cos_alpha, cos_beta, amp + + +class Sphere: + def __init__(self, nt): + """Sphere class. nt is the number of spherical triangles along the edge of the + sphere.""" + self.nt = nt + + def coordinates(self, nt: int): + """Coordinates over the position octant face of the sphere. + + Args: + nt: Number of spherical triangles along the edge of the spherical face. + """ + + # Equal area triangle mapping from octahedron to sphere. + # Reference: Computers and Mathematics with Applications 67 (2014) 1092–1107 + + octa = Octahedron(nt) + x, y, z = octa.coordinates(nt) + + z_adj = z / nt + factor1 = np.pi / 2 * (y / (x + y)) + factor1[-1] = 0 + z_sphere = 2 * z_adj * (1 - z_adj / 2) + + factor2 = np.sqrt(1 - z_sphere ** 2) + factor2[-1] = 0 + + x_sphere = factor2 * np.cos(factor1) + y_sphere = factor2 * np.sin(factor1) + + return x_sphere, y_sphere, z_sphere + + def cosines_and_amplitudes(self, integration_density: int = 72): + r""" + Calculates and return the direction cosines and the related amplitudes for + the positive quadrant of the sphere. The direction cosines corresponds to + angle $\alpha$ and $\beta$, where $\alpha$ is the azimuthal angle and + $\beta$ is the polar angle. The amplitudes are evaluated as $\frac{1}{r^3}$ + where $r$ is the distance from the origin to the face of the unit + octahedron in the positive quadrant along the line given by the values of + $\alpha$ and $\beta$. + + :ivar integration_density: + The value is an integer which represents the frequency of class I + geodesic polyhedra. These polyhedra are used in calculating the + spherical average. Presently we only use octahedral as the frequency1 + polyhedra. As the frequency of the geodesic polyhedron increases, the + polyhedra approach a sphere geometry. A higher frequency will result in + a better powder averaging. The default value is 72. + Read more on the `Geodesic polyhedron + `_. + + :return cos_alpha: The cosine of the azimuthal angle. + :return cos_beta: The cosine of the polar angle. + :return amp: The amplitude at the given $\alpha$ and $\beta$. + """ + nt = integration_density + x, y, z = self.coordinates(nt) + + cos_beta = z + cos_alpha = np.cos(0.5 * np.pi * (y / (x + y))) + amp = np.ones(cos_beta.size) + return cos_alpha, cos_beta, amp def triangle_interpolation(f, spec, amp=1.0): diff --git a/tests/test_lineshapes.py b/tests/test_lineshapes.py index 37062690b..3d69a0f6a 100644 --- a/tests/test_lineshapes.py +++ b/tests/test_lineshapes.py @@ -228,7 +228,7 @@ def test_pure_quadrupolar_lineshapes(): filename = path.join(path_, f"test{i:02d}", f"test{i:02d}.json") data_mrsimulator, data_source = c_setup(*get_data(filename)) np.testing.assert_almost_equal( - data_mrsimulator, data_source, decimal=2, err_msg=message + data_mrsimulator, data_source, decimal=1, err_msg=message ) # random euler angle. Euler angles should not affect the spectrum. diff --git a/tests/test_orientations_and_triangle_interpolation.py b/tests/test_orientations_and_triangle_interpolation.py index bfc9b4e92..e0bec5c20 100644 --- a/tests/test_orientations_and_triangle_interpolation.py +++ b/tests/test_orientations_and_triangle_interpolation.py @@ -3,14 +3,19 @@ import mrsimulator.tests.tests as clib import numpy as np -from .python_test_for_c_code.orientation import cosine_of_polar_angles_and_amplitudes +from .python_test_for_c_code.orientation import Octahedron from .python_test_for_c_code.orientation import triangle_interpolation +# import matplotlib.pyplot as plt + def test_octahedron_averaging_setup(): nt = 64 - cos_alpha_py, cos_beta_py, amp_py = cosine_of_polar_angles_and_amplitudes(nt) - exp_I_alpha_c, exp_I_beta_c, amp_c = clib.cosine_of_polar_angles_and_amplitudes(nt) + octa = Octahedron(nt) + cos_alpha_py, cos_beta_py, amp_py = octa.cosines_and_amplitudes(nt) + exp_I_alpha_c, exp_I_beta_c, amp_c = clib.cosine_of_polar_angles_and_amplitudes( + nt, integration_type=0 + ) assert np.allclose(cos_alpha_py, exp_I_alpha_c.real, atol=1e-15) assert np.allclose(cos_beta_py, exp_I_beta_c.real, atol=1e-15) @@ -38,6 +43,12 @@ def test_triangle_interpolation(): amp_c = np.zeros(100) clib.triangle_interpolation(list_, amp_c) + # plt.plot(amp_py, 'b', label='py') + # plt.plot(amp_c, 'r--', label='c') + # plt.plot(list_, [0, 2/(list_[2]-list_[0]), 0], 'k*', label='c') + # plt.legend() + # plt.show() + assert np.allclose(amp_py, amp_c, atol=1e-15) @@ -66,6 +77,13 @@ def test_triangle_rasterization(): clib.triangle_interpolation(lst1, amp2) clib.triangle_interpolation(lst2, amp3) + # plt.imshow(amp1, origin='lower', cmap='gray', aspect='auto') + # x_t, y_t = np.append(lst2, lst2[0]) - 0.5, np.append(lst1, lst1[0]) - 0.5 + # plt.plot(x_t, y_t, 'r', label='vertex') + # plt.scatter(lst2-0.5, lst1-0.5, s=20, color='yellow', label='vertex') + # plt.legend() + # plt.show() + assert np.allclose(amp2, amp1.sum(axis=1), atol=1e-15) assert np.allclose(amp3, amp1.sum(axis=0), atol=1e-15) @@ -85,6 +103,13 @@ def test_triangle_rasterization(): clib.triangle_interpolation2D(lst1, lst2, amp1) clib.triangle_interpolation(lst1, amp2) + # plt.imshow(amp1, origin='lower', cmap='gray', aspect='auto') + # x_t, y_t = np.append(lst2, lst2[0]) - 0.5, np.append(lst1, lst1[0]) - 0.5 + # plt.plot(x_t, y_t, 'r', label='vertex') + # plt.scatter(lst2-0.5, lst1-0.5, s=20, color='yellow', label='vertex') + # plt.legend() + # plt.show() + assert np.allclose(amp2, amp1.sum(axis=1), atol=1e-15) # triangles with one or more vertices outside a grid voxel @@ -101,4 +126,11 @@ def test_triangle_rasterization(): clib.triangle_interpolation2D(lst1, lst2, amp1) clib.triangle_interpolation(lst1, amp2) + # plt.imshow(amp1, origin='lower', cmap='gray', aspect='auto') + # x_t, y_t = np.append(lst2, lst2[0]) - 0.5, np.append(lst1, lst1[0]) - 0.5 + # plt.plot(x_t, y_t, 'r', label='vertex') + # plt.scatter(lst2-0.5, lst1-0.5, s=20, color='yellow', label='vertex') + # plt.legend() + # plt.show() + assert np.allclose(amp2, amp1.sum(axis=1), atol=1e-15)