From 47feebe22fc6a8bfad8dc3a2e0cbf285c31b67a9 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Fri, 2 May 2025 11:35:01 +0000 Subject: [PATCH] use C function and integration for anisotropic diffusion --- src/nmreval/clib/diffusion.c | 16 +++++++++++++++ src/nmreval/clib/diffusion.so | Bin 0 -> 15840 bytes src/nmreval/distributions/helper.py | 13 ++++++++++-- src/nmreval/models/diffusion.py | 30 +++++++++++++++++++++++----- 4 files changed, 52 insertions(+), 7 deletions(-) create mode 100644 src/nmreval/clib/diffusion.c create mode 100755 src/nmreval/clib/diffusion.so diff --git a/src/nmreval/clib/diffusion.c b/src/nmreval/clib/diffusion.c new file mode 100644 index 0000000..a1c71b6 --- /dev/null +++ b/src/nmreval/clib/diffusion.c @@ -0,0 +1,16 @@ +/* integrands used in quadrature integration with scipy's LowLevelCallables */ +#include + +double anistropicDiffusion(double x, void *user_data) { + double *c = (double *)user_data; + + double q = c[0]; + double t = c[1]; + double d_perp = c[2]; + double d_par = c[3]; + + double cos_theta = cos(x); + double sin_theta = sin(x); + + return exp(-q * q * t * (d_par * cos_theta * cos_theta + d_perp * sin_theta * sin_theta)) * sin_theta; +} \ No newline at end of file diff --git a/src/nmreval/clib/diffusion.so b/src/nmreval/clib/diffusion.so new file mode 100755 index 0000000000000000000000000000000000000000..8bc6c6fa08b2a0ca684c70a756b33e88c3b4571c GIT binary patch literal 15840 zcmeHOU2Ggz6~60@Q|HgSDGs;^C|i_>q>zoBDygALvW~sZw2D*LwMimT&5pgh_R{@v zcW0>`N^2k5SwB&<#;wico;k`PDgF)5Hrs3EnD zu6xx(PVAOVS4-)ZIH8r)%WZUXwAexYaF<&ZH2lnSs*JgQ70^Ey7~0F*O1_D>8!&Oo zuh+R^lDcT{xEst%2`?o)=w*qA?%_&ExR3l9lKTki0~|yG&9N*KAK7UxM`LWYYQ#f> z_cZaipZ~5P?UT}8r@XAl{Xru=+-@Iz!t^zX8(p=j3D)sH+Fo9btf4kN34yLN$v0);kY+1t*`&t1+at4QQ>X(2Z|n#*RZm0YQ)W=j=S$rV**aemFRDxZ|pOaOe3J!iEQb0F{vLPq+qSJeMG`gj-@9 zbDWVw3BD@Dcym5Eh`#6y=nUu#=nUu#=nUu#=nUu#=nVX;3>a5F?lG=)|M93&#<^(BlNEE{(ZZTD)^jK(9@`8Bl|};`O7DtCYP!noi?tV{DfcPjjJEoJB@26j}gB9 zEOl{xo-XfpA0|8WEZT*>^G$M26NeHtaWAqvdBeE+trg>1a>cls+#tPiFA`PksEd*4 zVQU>$l&u=9eztCtbt7v;iZ|DPdcV=2(n~j~wzW<8Ef0K)8c5dv9C>NiL@#ABa$)se zRBj5GIZ^tn!LI{!l*?u~w@Yh>5fu8*k&j?khrpfjK| zpfjK|pfjK|pfjK|pfjK|pfm8lkOAy7#lF&==#=a;-6eMagKps=v9VuupV)ZS+bcHq zt@8UJ3;ch%-zc%2l1-@CbDEOPq}U_+UGZ<{V3Mk<)DF@8qJ_qr+#L~`GQ^MfUhHB0 zjNr>#+ye01QW5(3V~N`)8twh9;r74TuK5333j5fD?VXVE8Wa6h(byM!c4Xu$u|uPo znVeOO4GzSQ#Rrd^I4q^uv4LZQ1A~JD!CL)xn06FnU%b~7rd{J0Pd^@37)L+8MPa=C zc(=kh`0=d@vG3a-|EcD4;m7+F=7k^MqlOY*#4U%_r&~bY zK(a;MZ3G?s{n17H>Twcfdeusd3ISe8$eC z;{)v@`0z%~w8~{`$;=e(@{-DytwP33R||zDqIfKmB5hx28mBy{Q>HmKIXs>;ljleI zILYXx^TXp$j1WvmNdTJZO39qJifKMn!nKi`4(XUlL&(PHq;h(s%(Lg7Ix~FEd}?g$ zndFo?HGJk=lB)HeH*pTFv^|)DQz!rISPP#>X?J2Ja56^4D@z62njvkMT|JLdv1DiB zbH!?WrkcyAkL1$IvF5GHyo#rnio|raT~_gOCT}q)rG>n$;!YCcRE*D+$h0$ybmbwT zpi_uGlakFfQwN-SW-1jGr}Uulg2eo6|^;A8!P^$IE&q5`jxz{l?|^lb^`cIlC7 zuFJ5#Nf5?$gOA^7Xsp91#1XE;E4=RFy1~cqIW*Q)!TzKDN#UQ7O7Xi7&GoY|pT5mn z`p^H1;67m;+5EDSFDKL%-RHT7_OV{W zdDLYYKioHX$YYkc90ESpb+w}elEj?r<^_D{8|3oxRYGD1TrEn#*8Z=O$IHihWI%{$ z2NeYEUy*|`^dIZVp$>d*%Pa4Z!C0`L;D5mzEU^rs9@dyUj6+gV;ff@lzi%F)5BiVq jAAUa2cW_xpkGs}oG$hV|Nh`kg{8vu6gb#%W73BW|LgSB( literal 0 HcmV?d00001 diff --git a/src/nmreval/distributions/helper.py b/src/nmreval/distributions/helper.py index b165f27..4927831 100644 --- a/src/nmreval/distributions/helper.py +++ b/src/nmreval/distributions/helper.py @@ -5,6 +5,17 @@ from ctypes import CDLL, c_double, c_void_p from ..lib.logger import logger +diffusion_lib = None +try: + diffusion_lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'diffusion.so')) + + diffusion_lib.anistropicDiffusion.restype = c_double + diffusion_lib.anistropicDiffusion.argtypes = (c_double, c_void_p) + + HAS_C_FUNCS = True +except OSError: + HAS_C_FUNCS = False + lib = None try: lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so')) @@ -39,10 +50,8 @@ try: lib.energyDistSuscImag.restype = c_double lib.energyDistSuscImag.argtypes = (c_double, c_void_p) - HAS_C_FUNCS = True logger.info('Use C functions') except OSError: HAS_C_FUNCS = False logger.info('Use python functions') - diff --git a/src/nmreval/models/diffusion.py b/src/nmreval/models/diffusion.py index b84d08e..fc492f7 100644 --- a/src/nmreval/models/diffusion.py +++ b/src/nmreval/models/diffusion.py @@ -1,7 +1,11 @@ +from ctypes import c_double, cast, c_void_p, pointer + import numpy as np -from scipy import special as special +from scipy import special as special, LowLevelCallable +from scipy.integrate import quad from ..utils import gamma +from nmreval.distributions.helper import HAS_C_FUNCS, diffusion_lib class Diffusion: @@ -103,13 +107,29 @@ class AnisotropicDiffusion(object): tp = x relax = np.exp(-(tp/trel)**brel)*np.exp(-(tp/trel)**brel) - q_squared = np.power(g * nucleus * tp, 2) + q = g * nucleus * tp t = 2 * tp / 3 + tm - z = np.sqrt(q_squared * (d_par - d_perp) * t) # Callaghan eq (6.89) - diffs = np.exp(-q_squared*t*d_perp) * special.erf(z) / z + if HAS_C_FUNCS: + # divide by 2 to normalize by integral sin(x), x=0..pi + diffusion_decay = AnisotropicDiffusion._integrate_c(q, t, d_perp, d_par) / 2 + else: + z = np.sqrt(q**2 * (d_par - d_perp) * t) + diffusion_decay = np.exp(-q**2 * t * d_perp) * special.erf(z) / z - return m0 * diffs * relax + return m0 * diffusion_decay * relax + + @staticmethod + def _integrate_c(q, t, d_perp, d_par) -> np.ndarray: + diffusion_decay = np.zeros_like(t) + + for (i, t_i) in enumerate(t): + c = (c_double * 4)(q, t_i, d_perp, d_par) + user_data = cast(pointer(c), c_void_p) + + diffusion_decay[i] = quad(LowLevelCallable(diffusion_lib.anistropicDiffusion, user_data), 0, np.pi, epsabs=1e-13)[0] + + return diffusion_decay class Peschier: