IPKM/nmreval
16
2
Fork 1
Code Issues 41 Pull Requests Actions Packages Projects Releases Wiki Activity

add C function to integrate anisotropic diffusion

Browse Source
This commit is contained in:
Dominik Demuth 2025-05-02 13:25:56 +02:00
parent 81357a3e4f
commit e4c54b7104
2 changed files with 13 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/nmreval/distributions/helper.py
View File

@ -55,4 +55,3 @@ try:
except OSError:
HAS_C_FUNCS = False
logger.info('Use python functions')

35
src/nmreval/models/diffusion.py
View File

@ -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: