changing default order of hdf works with node names and displays correct label (#310)

Fix problem with incorrect selection of label and group of HDF files when using not default order

Co-authored-by: Dominik Demuth <dominik.demuth@physik.tu-darmstadt.de>
Reviewed-on: #310

src/nmreval/clib/diffusion.c

@@ -0,0 +1,16 @@
+/* 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,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')
-

src/nmreval/io/hdfreader.py

@@ -302,15 +302,26 @@ class HdfReader(HdfNode):
 
     def make_signal(self, node, flag: str = 'fid', value: str = None, group: str = None):
         if value is None:
+            data_name = node.name
             value = self._get_parameter_values(node, node.parameter)
         else:
             try:
+                data_name = f"{value}={node.parameter[value]}"
                 value = node.parameter[value]
             except KeyError:
+                print(node.title_parameter)
                 try:
-                    value = node.title_parameter[1][value]
+                    temp = node
+                    while value != temp.title_parameter[0][0]:
+                        if temp.parent is None:
+                            break
+                        temp = temp.parent
+
+                    value = temp.title_parameter[0][1]
+                    data_name = temp.name
                 except KeyError:
                     print(f'{value} is not a valid key for {node.name}')
+                    data_name = node.name
                     value = None
 
         if group is None:
@@ -343,11 +354,11 @@ class HdfReader(HdfNode):
         dw = float(index['dwelltime'])
         if flag == 'fid':
             x = np.arange(len(y)) * dw
-            ret = FID(x, y, name=node.name, value=value, group=group, filename=self.file.filename)
+            ret = FID(x, y, name=data_name, value=value, group=group, filename=self.file.filename)
 
         elif flag == 'spectrum':
             x = np.linspace(-1/dw, 1/dw, num=len(y))
-            ret = Spectrum(x, y, name=node.name, value=value, group=group, filename=self.file.filename)
+            ret = Spectrum(x, y, name=data_name, value=value, group=group, filename=self.file.filename)
         else:
             raise ValueError(f'{flag} unknown, use `fid` or `spectrum`.')
 

src/nmreval/models/basic.py

@@ -203,6 +203,31 @@ class Sinc:
         return c * np.sinc(((x-x0)/w)/np.pi)
 
 
+class Sigmoid:
+    type = 'Basic'
+    name = 'Sigmoid'
+    equation = 'C / [1 + exp(-a * (x - x_{0})] + y_{0}'
+    params = ['C', 'a', 'x_{0}', 'y_{0}']
+
+    @staticmethod
+    def func(x, c, a, x0, y0):
+        """
+        Sigmoid function
+
+        .. math::
+            y = C / [1 + exp(-a * (x - x_0))] + y_0
+
+        Args:
+            x (array_like): Input values
+            c (float): Prefactor
+            a (float): Steepness of the sigmoid
+            x0 (float): x position of the sigmoid's midpoint
+            y0 (float): y position of the sigmoid's midpoint
+
+        """
+        return c / (1 + np.exp(-a * (x - x0))) + y0
+
+
 class Sine:
     """
     Wavy sine function

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: