From 81357a3e4f02e02719d81e8b8fa9a85ea5d3518e Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Fri, 2 May 2025 13:20:35 +0200 Subject: [PATCH 1/3] add C function to integrate anisotropic diffusion --- src/nmreval/clib/diffusion.c | 16 +++++++++++++++ src/nmreval/clib/diffusion.so | Bin 0 -> 15840 bytes src/nmreval/distributions/helper.py | 12 ++++++++++- src/nmreval/models/diffusion.py | 30 +++++++++++++++++++++++++++- 4 files changed, 56 insertions(+), 2 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..1138cdd 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,7 +50,6 @@ 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: diff --git a/src/nmreval/models/diffusion.py b/src/nmreval/models/diffusion.py index b84d08e..f43594f 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: @@ -111,6 +115,30 @@ class AnisotropicDiffusion(object): return m0 * diffs * relax + @staticmethod + def _integrate_c(x, m0, d_perp, d_par, trel, brel, g, t2, nucleus=2.67522128e8, t_axis='tm') -> np.ndarray: + diffusion_decay = np.zeros_like(x) + + if t_axis == 'tm': + tm = x + tp = t2 + relax = np.exp(-(tm/trel)**brel) + else: + tm = t2 + tp = x + relax = np.exp(-(tp/trel)**brel)*np.exp(-(tp/trel)**brel) + + q = g * nucleus * tp + t = 2 * tp / 3 + tm + + 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 m0 * diffusion_decay * relax + class Peschier: name = 'Diffusion + Cross-Relaxation' -- 2.39.5 From e4c54b7104cf8181798a6a4ea01cd69fdabf7b23 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Fri, 2 May 2025 13:25:56 +0200 Subject: [PATCH 2/3] add C function to integrate anisotropic diffusion --- src/nmreval/distributions/helper.py | 1 - src/nmreval/models/diffusion.py | 35 +++++++++++------------------ 2 files changed, 13 insertions(+), 23 deletions(-) diff --git a/src/nmreval/distributions/helper.py b/src/nmreval/distributions/helper.py index 1138cdd..4927831 100644 --- a/src/nmreval/distributions/helper.py +++ b/src/nmreval/distributions/helper.py @@ -55,4 +55,3 @@ try: 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 f43594f..c479d06 100644 --- a/src/nmreval/models/diffusion.py +++ b/src/nmreval/models/diffusion.py @@ -98,27 +98,6 @@ class AnisotropicDiffusion(object): @staticmethod def func(x, m0, d_perp, d_par, trel, brel, g, t2, nucleus=2.67522128e8, t_axis='tm'): - if t_axis == 'tm': - tm = x - tp = t2 - relax = np.exp(-(tm/trel)**brel) - else: - tm = t2 - tp = x - relax = np.exp(-(tp/trel)**brel)*np.exp(-(tp/trel)**brel) - - q_squared = np.power(g * nucleus * tp, 2) - 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 - - return m0 * diffs * relax - - @staticmethod - def _integrate_c(x, m0, d_perp, d_par, trel, brel, g, t2, nucleus=2.67522128e8, t_axis='tm') -> np.ndarray: - diffusion_decay = np.zeros_like(x) - if t_axis == 'tm': tm = x tp = t2 @@ -130,6 +109,18 @@ class AnisotropicDiffusion(object): q = g * nucleus * tp t = 2 * tp / 3 + tm + # Callaghan eq (6.89) + if HAS_C_FUNCS: + diffusion_decay = AnisotropicDiffusion._integrate_c(q, t, d_perp, d_par) + 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 * 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) @@ -137,7 +128,7 @@ class AnisotropicDiffusion(object): diffusion_decay[i] = quad(LowLevelCallable(diffusion_lib.anistropicDiffusion, user_data), 0, np.pi, epsabs=1e-13)[0] - return m0 * diffusion_decay * relax + return diffusion_decay class Peschier: -- 2.39.5 From e5a74f3af2c34d740117396c9a405560f81b48f9 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Fri, 2 May 2025 11:32:51 +0000 Subject: [PATCH 3/3] Update src/nmreval/models/diffusion.py --- src/nmreval/models/diffusion.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nmreval/models/diffusion.py b/src/nmreval/models/diffusion.py index c479d06..fc492f7 100644 --- a/src/nmreval/models/diffusion.py +++ b/src/nmreval/models/diffusion.py @@ -111,7 +111,8 @@ class AnisotropicDiffusion(object): t = 2 * tp / 3 + tm # Callaghan eq (6.89) if HAS_C_FUNCS: - diffusion_decay = AnisotropicDiffusion._integrate_c(q, t, d_perp, d_par) + # 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 -- 2.39.5