add sigmoid function to basic models

src/nmreval/clib/diffusion.c

@@ -1,16 +0,0 @@
-/* integrands used in quadrature integration with scipy's LowLevelCallables */
-#include <math.h>
-
-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;
-}

src/nmreval/distributions/helper.py

@@ -5,17 +5,6 @@ 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'))
@@ -50,8 +39,10 @@ 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')
+

src/nmreval/models/diffusion.py

@@ -1,11 +1,7 @@
-from ctypes import c_double, cast, c_void_p, pointer
-
 import numpy as np
-from scipy import special as special, LowLevelCallable
+from scipy import special as special
-from scipy.integrate import quad
 
 from ..utils import gamma
-from nmreval.distributions.helper import HAS_C_FUNCS, diffusion_lib
 
 
 class Diffusion:
@@ -107,29 +103,13 @@ class AnisotropicDiffusion(object):
             tp = x
             relax = np.exp(-(tp/trel)**brel)*np.exp(-(tp/trel)**brel)
 
-        q = g * nucleus * tp
+        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)
-        if HAS_C_FUNCS:
+        diffs = np.exp(-q_squared*t*d_perp) * special.erf(z) / z
-            # 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 * diffusion_decay * relax
+        return m0 * diffs * 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: