merge cfunc -> main

AppImageBuilder.yml

@@ -36,10 +36,10 @@ AppDir:
 
     include:
       # for /usr/bin/env
-      - coreutils
-      - dash
-      - zsync
-      - hicolor-icon-theme
+#      - coreutils
+#      - dash
+#      - zsync
+#      - hicolor-icon-theme
       - libatlas3-base
       - python3.9-minimal
       - python3-numpy
@@ -57,17 +57,24 @@ AppDir:
       - libqt5test5
       - libqt5xml5
       - qtbase5-dev-tools
-      - qtchooser
       - pyqt5-dev-tools
       - libavahi-client3
       - libavahi-common-data
       - libavahi-common3
       - libwacom2
       - libwacom-common
-  after_bundle: |
-    echo "MONSTER SED FOLLOWING...(uncomment if needed for mpl-data)"
-    # sed -i s,\'/usr/share/matplotlib/mpl-data\',"f\"\{os.environ.get\('APPDIR'\,'/'\)\}/usr/share/matplotlib/mpl-data\"", ${TARGET_APPDIR}/usr/lib/python3/dist-packages/matplotlib/__init__.py
+
+  files:
+    exclude:
+      - usr/share/man
+      - usr/share/doc/*/README.*
+      - usr/share/doc/*/changelog.*
+      - usr/share/doc/*/NEWS.*
+      - usr/share/doc/*/TODO.}*
   runtime:
+    # if needed, apparently replaces hardcoded location with APPDIR location
+    # path_mappings:
+    #   - /usr/share/matplotlib/mpl-data:$APPDIR/usr/share/matplotlib/mpl-data
     version: "continuous"
     env:
       PATH: '${APPDIR}/usr/bin:${PATH}'

Makefile

@@ -8,28 +8,34 @@ PYRCC = pyrcc5
 RESOURCE_DIR = src/resources/_ui
 
 #Directory for compiled resources
-COMPILED_DIR = src/gui_qt/_py
+PYQT_DIR = src/gui_qt/_py
 
 #UI files to compile, uses every *.ui found in RESOURCE_DIR
 UI_FILES = $(foreach dir, $(RESOURCE_DIR), $(notdir $(wildcard $(dir)/*.ui)))
-COMPILED_UI = $(UI_FILES:%.ui=$(COMPILED_DIR)/%.py)
+PYQT_UI = $(UI_FILES:%.ui=$(PYQT_DIR)/%.py)
 
 SVG_FILES = $(foreach dir, $(RCC_DIR), $(notdir $(wildcard $(dir)/*.svg)))
 PNG_FILES = $(SVG_FILES:%.svg=$(RCC_DIR)/%.png)
 
-all : ui
- 
-ui : $(COMPILED_UI)
+CC = gcc
+CFLAGS = -O2 -fPIC
+LDFLAGS = -shared
 
-rcc:  $(PNG_FILES)
+C_DIR = src/nmreval/clib
+
+all : ui compile
+
+ui : $(PYQT_UI)
+
+rcc : $(PNG_FILES)
+
+# only one C file at the moment
+compile : $(C_DIR)/integrate.c
+	$(CC) $(LDFLAGS) $(CFLAGS) -o $(C_DIR)/integrate.so $<
 
- 
 $(COMPILED_DIR)/%.py : $(RESOURCE_DIR)/%.ui
 	$(PYUIC) $< -o $@
-# 	replace import of ressource to correct path
-#	@sed -i s/images_rc/nmrevalqt.$(COMPILED_DIR).images_rc/g $@
-#	@sed -i /images_rc/d $@
-	
+
 $(RCC_DIR)/%.png : $(RCC_DIR)/%.svg
 	convert -background none $< $@
 	$(PYRCC) $(RCC_DIR)/images.qrc -o $(COMPILED_DIR)/images_rc.py

src/nmreval/clib/integrate.c

@@ -0,0 +1,158 @@
+/* integrands used in quadrature integration with scipy's LowLevelCallables */
+#include <math.h>
+
+const double KB = 8.617333262145179e-05;
+
+/* FFHS functions */
+double ffhsSD(double x, void *user_data) {
+	double *c = (double *)user_data;
+
+	double omega = c[0];
+	double tau = c[1];
+
+	double u = x*x;
+	double res = u*u * tau;
+
+	res /= 81. + 9.*u - 2.*u*u + u*u*u;
+	res /= u*u + omega*omega * tau*tau;
+	return res;
+}
+
+/* log-gaussian functions */
+double logNormalDist(double tau, double tau0, double sigma) {
+    return exp(- pow((log(tau/tau0) / sigma), 2) / 2.) / sqrt(2*M_PI)/sigma;
+}
+
+double logGaussian_imag_high(double u, void *user_data) {
+	double *c = (double *)user_data;
+
+	double omega = c[0];
+	double tau = c[1];
+	double sigma = c[2];
+
+    double uu = exp(-u);
+    double dist = logNormalDist(1./uu, tau, sigma);
+
+    return dist * omega * uu / (pow(uu, 2) + pow(omega, 2));
+}
+
+double logGaussian_imag_low(double u, void *user_data) {
+	double *c = (double *)user_data;
+
+	double omega = c[0];
+	double tau = c[1];
+	double sigma = c[2];
+
+    double uu = exp(u);
+
+    double dist = logNormalDist(uu, tau, sigma);
+
+    return dist * omega * uu / (1. + pow(omega*uu, 2));
+}
+
+double logGaussian_real_high(double u, void *user_data) {
+	double *c = (double *)user_data;
+
+	double omega = c[0];
+	double tau = c[1];
+	double sigma = c[2];
+
+    double uu = exp(-2.*u);
+    double dist = logNormalDist(exp(uu), tau, sigma);
+
+    return dist * uu / (uu + pow(omega, 2));
+}
+
+double logGaussian_real_low(double u, void *user_data) {
+	double *c = (double *)user_data;
+
+	double omega = c[0];
+	double tau = c[1];
+	double sigma = c[2];
+
+    double uu = exp(u);
+
+    double dist = logNormalDist(uu, tau, sigma);
+
+    return dist / (1. + pow(omega*uu, 2));
+}
+
+double logGaussianCorrelation(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double t = c[0];
+	double tau = c[1];
+	double sigma = c[2];
+
+    double uu = exp(x);
+
+    double dist = logNormalDist(uu, tau, sigma);
+
+    return dist * exp(-t/uu);
+}
+
+// functions for distribution of energy
+double normalDist(double x, double x0, double sigma) {
+    return exp(- pow((x-x0) / sigma, 2) / 2.) / sqrt(2 * M_PI) / sigma;
+}
+
+double rate(double tau0, double ea, double t) {
+    return exp(-ea / t / KB) / tau0;
+}
+
+double energyDist_SD(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double omega = c[0];
+    double tau0 = c[1];
+    double e_m = c[2];
+    double e_b = c[3];
+    double temp = c[4];
+
+    double r = rate(tau0, x, temp);
+
+    return r/(pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
+}
+
+double energyDistSuscReal(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double omega = c[0];
+    double tau0 = c[1];
+    double e_m = c[2];
+    double e_b = c[3];
+    double temp = c[4];
+
+    double r = rate(tau0, x, temp);
+
+    return 1 / (pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
+}
+
+double energyDistSuscImag(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double omega = c[0];
+    double tau0 = c[1];
+    double e_m = c[2];
+    double e_b = c[3];
+    double temp = c[4];
+
+    double r = rate(tau0, x, temp);
+
+    return omega * r / (pow(r, 2) + pow(omega, 2)) * normalDist(x, e_m, e_b);
+}
+
+double energyDistCorrelation(double x, void *user_data) {
+    double *c = (double *)user_data;
+
+    double t = c[0];
+    double tau0 = c[1];
+    double e_m = c[2];
+    double e_b = c[3];
+    double temp = c[4];
+
+    double r = rate(tau0, x, temp);
+
+    return normalDist(x, e_m, e_b) * exp(-t * r);
+}
+

src/nmreval/distributions/base.py

@@ -29,7 +29,7 @@ class Distribution(abc.ABC):
 
     @staticmethod
     @abc.abstractmethod
-    def susceptibility(omega, tau, *args):
+    def susceptibility(omega: ArrayLike, tau: ArrayLike, *args: Any):
         pass
 
     @classmethod

src/nmreval/distributions/energy.py

@@ -1,13 +1,18 @@
 from itertools import product
+from ctypes import c_double, cast, pointer, c_void_p
 
 import numpy as np
+from scipy import LowLevelCallable
 from scipy.integrate import quad, simps as simpson
 
 from .base import Distribution
 from ..lib.utils import ArrayLike
 from ..utils.constants import kB
 
+from .helper import HAS_C_FUNCS, lib
 
+
+# noinspection PyMethodOverriding
 class EnergyBarriers(Distribution):
     name = 'Energy barriers'
     parameter = [r'\tau_{0}', r'E_{m}', r'\Delta E']
@@ -23,24 +28,30 @@ class EnergyBarriers(Distribution):
         return np.exp(-e_a / (kB * te)) / t0
 
     @staticmethod
-    def energydistribution(e_a, mu, sigma):
+    def energy_distribution(e_a, mu, sigma):
         return np.exp(-0.5 * ((mu-e_a) / sigma) ** 2) / (np.sqrt(2 * np.pi) * sigma)
 
     @staticmethod
-    def correlation(t, temperature, *args):
-        tau0, e_m, e_b = args
-
-        def integrand(e_a, ti, t0, mu, sigma, te):
-            # correlation time would go to inf for higher energies, so we use rate
-            return np.exp(-ti*EnergyBarriers.rate(t0, e_a, te)) * EnergyBarriers.energydistribution(e_a, mu, sigma)
-
+    def correlation(t: ArrayLike, temperature: ArrayLike, tau0: float, e_m: float, e_b: float) -> ArrayLike:
         t = np.atleast_1d(t)
         temperature = np.atleast_1d(temperature)
 
-        e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
+        if HAS_C_FUNCS:
+            ret_val = _integrate_c(lib.energyDistCorrelation, t, temperature, tau0, e_m, e_b)
+        else:
+            ret_val = _integrate_py(_integrand_time, t, temperature, tau0, e_m, e_b)
 
-        ret_val = np.array([simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis)
-                            for o in t for tt in temperature])
+        return ret_val
+
+    @staticmethod
+    def specdens(omega: ArrayLike, temperature: ArrayLike, tau0: float, e_m: float, e_b: float) -> ArrayLike:
+        omega = np.atleast_1d(omega)
+        temperature = np.atleast_1d(temperature)
+
+        if HAS_C_FUNCS:
+            ret_val = _integrate_c(lib.energyDist_SD, omega, temperature, tau0, e_m, e_b)
+        else:
+            ret_val = _integrate_py(_integrand_sd, omega, temperature, tau0, e_m, e_b)
 
         return ret_val
 
@@ -51,70 +62,73 @@ class EnergyBarriers(Distribution):
         omega = np.atleast_1d(omega)
         temperature = np.atleast_1d(temperature)
 
-        e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
-        ret_val = []
-        for o, tt in product(omega, temperature):
-            ret_val.append(simpson(_integrand_freq_real(e_axis, o, tau0, e_m, e_b, tt), e_axis) -
-                           1j * simpson(_integrand_freq_imag(e_axis, o, tau0, e_m, e_b, tt), e_axis))
-
-        return np.array(ret_val)
-
-    @staticmethod
-    def specdens(omega, temperature, *args):
-        # in contrast to other spectral densities, it's omega and temperature
-        tau0, e_m, e_b = args
-
-        def integrand(e_a, w, t0, mu, sigma, t):
-            r = EnergyBarriers.rate(t0, e_a, t)
-            return r/(r**2 + w**2) * EnergyBarriers.energydistribution(e_a, mu, sigma)
-
-        omega = np.atleast_1d(omega)
-        temperature = np.atleast_1d(temperature)
-
-        e_axis = np.linspace(max(0, e_m-50*e_b), e_m+50*e_b, num=5001)
-
-        ret_val = np.array([simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis)
-                            for o in omega for tt in temperature])
+        if HAS_C_FUNCS:
+            ret_val = _integrate_c(lib.energyDistSuscReal, omega, temperature, tau0, e_m, e_b) + \
+                      1j * _integrate_c(lib.energyDistSuscImag, omega, temperature, tau0, e_m, e_b)
+        else:
+            ret_val = _integrate_py(_integrand_susc_real, omega, temperature, tau0, e_m, e_b) + \
+                      1j * _integrate_py(_integrand_susc_imag, omega, temperature, tau0, e_m, e_b)
 
         return ret_val
 
     @staticmethod
-    def mean(*args):
-        return args[1]*np.exp(args[2]/(kB*args[0]))
+    def mean(temperature, tau0, ea):
+        return tau0*np.exp(ea/(kB*temperature))
 
     @staticmethod
-    def logmean(*args):
-        return args[1] + args[2] / (kB * args[0])
+    def logmean(temperature, tau0, ea):
+        return tau0 + ea / (kB * temperature)
 
     @staticmethod
     def max(*args):
         return args[1] * np.exp(args[2] / (kB * args[0]))
 
 
-def _integrate_process_imag(args):
-    pass
+# helper functions
+def _integrate_c(func, omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray:
+    res = []
+    for o, t in product(omega, temperature):
+        c = (c_double * 5)(o, tau0, e_m, e_b, t)
+        user_data = cast(pointer(c), c_void_p)
+        area = quad(LowLevelCallable(func, user_data), 0, np.infty, epsabs=1e-13)[0]
+
+        res.append(area)
+
+    ret_val = np.array(res).reshape(omega.shape[0], temperature.shape[0])
+
+    return ret_val.squeeze()
 
 
-def _integrate_process_real(args):
-    omega_i, t, tau0, mu, sigma, temp_j = args
-    return quad(_integrand_freq_real(), 0, 10, args=(omega_i, t, tau0, mu, sigma, temp_j))[0]
+def _integrate_py(func, axis, temp, tau0, e_m, e_b):
+    x = np.atleast_1d(axis)
+    temperature = np.atleast_1d(temp)
+
+    e_axis = np.linspace(max(0., e_m - 50*e_b), e_m + 50*e_b, num=5001)
+    ret_val = []
+    for o, tt in product(x, temperature):
+        ret_val.append(simpson(func(e_axis, o, tau0, e_m, e_b, tt), e_axis))
+
+    ret_val = np.array(ret_val).reshape(x.shape[0], temperature.shape[0])
+
+    return ret_val.squeeze()
 
 
-def _integrate_process_time(args):
-    omega_i, t, tau0, mu, sigma, temp_j = args
-    return quad(_integrand_time, 0, 10, args=(omega_i, t, tau0, mu, sigma, temp_j))[0]
-
-
-def _integrand_freq_real(u, omega, tau0, mu, sigma, temp):
+# python integrands
+def _integrand_sd(u, omega, tau0, mu, sigma, temp):
     r = EnergyBarriers.rate(tau0, u, temp)
-    return 1 / (r**2 + omega**2) * EnergyBarriers.energydistribution(u, mu, sigma)
+    return r / (r**2 + omega**2) * EnergyBarriers.energy_distribution(u, mu, sigma)
 
 
-def _integrand_freq_imag(u, omega, tau0, mu, sigma, temp):
+def _integrand_susc_real(u, omega, tau0, mu, sigma, temp):
+    r = EnergyBarriers.rate(tau0, u, temp)
+    return 1 / (r**2 + omega**2) * EnergyBarriers.energy_distribution(u, mu, sigma)
+
+
+def _integrand_susc_imag(u, omega, tau0, mu, sigma, temp):
     rate = EnergyBarriers.rate(tau0, u, temp)
-    return omega * rate / (rate**2 + omega**2) * EnergyBarriers.energydistribution(u, mu, sigma)
+    return omega * rate / (rate**2 + omega**2) * EnergyBarriers.energy_distribution(u, mu, sigma)
 
 
 def _integrand_time(u, t, tau0, mu, sigma, temp):
     rate = EnergyBarriers.rate(tau0, u, temp)
-    return EnergyBarriers.energydistribution(u, mu, sigma) * np.exp(-t*rate)
+    return EnergyBarriers.energy_distribution(u, mu, sigma) * np.exp(-t*rate)

src/nmreval/distributions/helper.py

@@ -0,0 +1,48 @@
+
+from pathlib import Path
+from ctypes import CDLL, c_double, c_void_p
+
+from ..lib.logger import logger
+
+
+lib = None
+try:
+    lib = CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so'))
+
+    # FFHS integrand
+    lib.ffhsSD.restype = c_double
+    lib.ffhsSD.argtypes = (c_double, c_void_p)
+
+    # Log-Gaussian integrands
+    lib.logGaussian_imag_high.restype = c_double
+    lib.logGaussian_imag_high.argtypes = (c_double, c_void_p)
+    lib.logGaussian_imag_low.restype = c_double
+    lib.logGaussian_imag_low.argtypes = (c_double, c_void_p)
+
+    lib.logGaussian_real_high.restype = c_double
+    lib.logGaussian_real_high.argtypes = (c_double, c_void_p)
+    lib.logGaussian_real_low.restype = c_double
+    lib.logGaussian_real_low.argtypes = (c_double, c_void_p)
+
+    lib.logGaussianCorrelation.restype = c_double
+    lib.logGaussianCorrelation.argtypes = (c_double, c_void_p)
+
+    # integrands for distribution of energies
+    lib.energyDist_SD.restype = c_double
+    lib.energyDist_SD.argtypes = (c_double, c_void_p)
+
+    lib.energyDistCorrelation.restype = c_double
+    lib.energyDistCorrelation.argtypes = (c_double, c_void_p)
+
+    lib.energyDistSuscReal.restype = c_double
+    lib.energyDistSuscReal.argtypes = (c_double, c_void_p)
+    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/distributions/intermolecular.py

@@ -1,9 +1,16 @@
+import ctypes
+
 import numpy as np
+from scipy import LowLevelCallable
 from scipy.integrate import quad
 
+from .helper import HAS_C_FUNCS, lib
 from .base import Distribution
 
 
+# Everything except spectral density is implemented in Python only because the only use case of FFHS is NMR
+# field cycling measurements with T1 results
+
 class FFHS(Distribution):
     name = 'Intermolecular (FFHS)'
     parameter = None
@@ -24,7 +31,7 @@ class FFHS(Distribution):
         return ret_val
 
     @staticmethod
-    def specdens(omega, tau0, *args):
+    def specdens_py(omega, tau0):
         def integrand(u, o, tau0):
             return u**4 * tau0 / (81 + 9*u**2 - 2*u**4 + u**6) / (u**4 + (o*tau0)**2)
             # return FFHS.distribution(u, tau0) * u / (1+o**2 * u**2)
@@ -33,6 +40,17 @@ class FFHS(Distribution):
 
         return ret_val * 54 / np.pi
 
+    @staticmethod
+    def specdens_c(omega, tau0):
+        res = []
+        for o in omega:
+            c = (ctypes.c_double * 2)(o, tau0)
+            user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
+            func = LowLevelCallable(lib.ffhsSD, user_data)
+            res.append(quad(func, 0, np.infty)[0])
+
+        return np.array(res) * 54 / np.pi
+
     @staticmethod
     def susceptibility(omega, tau0, *args):
         def integrand_real(u, o, tau0):
@@ -48,6 +66,9 @@ class FFHS(Distribution):
         return ret_val
 
 
+FFHS.specdens = FFHS.specdens_c if HAS_C_FUNCS else FFHS.specdens_py
+
+
 # class Bessel(Distribution):
 #     name = 'Intermolecular (Bessel)'
 #     parameter = None

src/nmreval/distributions/loggaussian.py

@@ -1,7 +1,12 @@
+import ctypes
 from multiprocessing import Pool, cpu_count
 from itertools import product
+from typing import Callable
 
 import numpy as np
+from scipy import LowLevelCallable
+
+from nmreval.lib.utils import ArrayLike
 
 try:
     from scipy.integrate import simpson
@@ -9,19 +14,21 @@ except ImportError:
     from scipy.integrate import simps as simpson
 from scipy.integrate import quad
 
-from .base import Distribution
+from nmreval.distributions.helper import HAS_C_FUNCS, lib
+from nmreval.distributions.base import Distribution
 
 
 __all__ = ['LogGaussian']
 
 
+# noinspection PyMethodOverriding
 class LogGaussian(Distribution):
     name = 'Log-Gaussian'
     parameter = [r'\sigma']
     bounds = [(0, 10)]
 
     @staticmethod
-    def distribution(tau, tau0, sigma: float):
+    def distribution(tau: ArrayLike, tau0: ArrayLike, sigma: float) -> ArrayLike:
         return np.exp(-0.5*(np.log(tau/tau0)/sigma)**2)/np.sqrt(2*np.pi)/sigma
 
     @staticmethod
@@ -29,54 +36,87 @@ class LogGaussian(Distribution):
         _t = np.atleast_1d(t)
         _tau = np.atleast_1d(tau0)
 
-        pool = Pool(processes=min(cpu_count(), 4))
-        integration_ranges = [(omega_i, tau_j, sigma) for (omega_i, tau_j) in product(_t, _tau)]
+        if HAS_C_FUNCS:
+            res = _integrate_correlation_c(_t, _tau, sigma)
+        else:
+            res = _integration_parallel(_t, _tau, sigma, _integrate_process_time)
 
-        with np.errstate(divide='ignore'):
-            res = np.array(pool.map(_integrate_process_time, integration_ranges))
-        ret_val = res.reshape((_t.shape[0], _tau.shape[0]))
-
-        return ret_val.squeeze()
+        return res.squeeze()
 
     @staticmethod
     def susceptibility(omega, tau0, sigma: float):
         _omega = np.atleast_1d(omega)
         _tau = np.atleast_1d(tau0)
 
-        pool = Pool(processes=min(cpu_count(), 4))
-        integration_ranges = [(omega_i, tau_j, sigma) for (omega_i, tau_j) in product(_omega, _tau)]
+        if HAS_C_FUNCS:
+            res_real = _integrate_susc_real_c(_omega, _tau, sigma)
+            res_imag = _integrate_susc_imag_c(_omega, _tau, sigma)
+        else:
+            res_real = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag)
+            res_imag = _integration_parallel(_omega, _tau, sigma, _integrate_process_real)
 
-        with np.errstate(divide='ignore'):
-            res_real = np.array(pool.map(_integrate_process_imag, integration_ranges))
-            res_imag = np.array(pool.map(_integrate_process_real, integration_ranges))
-        ret_val = (res_real+1j*res_imag).reshape((_omega.shape[0], _tau.shape[0]))
-
-        return ret_val.squeeze()
+        return (res_real + 1j * res_imag).squeeze()
 
     @staticmethod
-    def specdens(omega, tau0, sigma):
+    def specdens(omega: ArrayLike, tau: ArrayLike, sigma: float) -> np.ndarray:
         _omega = np.atleast_1d(omega)
-        _tau = np.atleast_1d(tau0)
+        _tau = np.atleast_1d(tau)
 
-        pool = Pool(processes=min(cpu_count(), 4))
-        integration_ranges = [(omega_i, tau_j, sigma) for (omega_i, tau_j) in product(_omega, _tau)]
-
-        with np.errstate(divide='ignore'):
-            res = np.array(pool.map(_integrate_process_imag, integration_ranges))
-        ret_val = res.reshape((_omega.shape[0], _tau.shape[0]))
+        if HAS_C_FUNCS:
+            ret_val = _integrate_susc_imag_c(_omega, _tau, sigma)
+        else:
+            ret_val = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag)
 
         ret_val /= _omega[:, None]
 
         return ret_val.squeeze()
 
-    def mean(*args):
-        return args[0]*np.exp(args[1]**2 / 2)
+    @staticmethod
+    def mean(tau, sigma):
+        return tau*np.exp(sigma**2 / 2)
+
+
+def _integration_parallel(x: np.ndarray, tau: np.ndarray, sigma: float, func: Callable) -> np.ndarray:
+    pool = Pool(processes=min(cpu_count(), 4))
+    integration_ranges = [(x_i, tau_j, sigma) for (x_i, tau_j) in product(x, tau)]
+
+    with np.errstate(divide='ignore'):
+        res = pool.map(func, integration_ranges)
+
+    res = np.array(res).reshape((x.shape[0], tau.shape[0]))
+
+    return res
+
+
+def _integrate_susc_imag_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> np.ndarray:
+    return _integrate_susc_c(lib.logGaussian_imag_low, lib.logGaussian_imag_high, omega, tau, sigma)
+
+
+def _integrate_susc_real_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> np.ndarray:
+    return _integrate_susc_c(lib.logGaussian_real_low, lib.logGaussian_real_high, omega, tau, sigma)
+
+
+def _integrate_susc_c(lowfunc, highfunc, omega, tau, sigma):
+    res = []
+
+    for o, t in product(omega, tau):
+        c = (ctypes.c_double * 3)(o, t, sigma)
+        user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
+
+        area = quad(LowLevelCallable(highfunc, user_data), 0, np.infty, epsabs=1e-12, epsrel=1e-12)[0]
+        area += quad(LowLevelCallable(lowfunc, user_data), -np.infty, 0, epsabs=1e-12, epsrel=1e-12)[0]
+
+        res.append(area)
+
+    res = np.asanyarray(res).reshape((omega.shape[0], tau.shape[0]))
+
+    return res
 
 
 def _integrate_process_imag(args):
     omega_i, tau_j, sigma = args
-    area = quad(_integrand_freq_imag_high, 0, 50, args=(omega_i, tau_j, sigma))[0]
-    area += quad(_integrand_freq_imag_low, -50, 0, args=(omega_i, tau_j, sigma))[0]
+    area = quad(_integrand_freq_imag_high, 0, 50, args=(omega_i, tau_j, sigma), epsabs=1e-12, epsrel=1e-12)[0]
+    area += quad(_integrand_freq_imag_low, -50, 0, args=(omega_i, tau_j, sigma), epsabs=1e-12, epsrel=1e-12)[0]
 
     return area
 
@@ -89,9 +129,25 @@ def _integrate_process_real(args):
     return area
 
 
+def _integrate_correlation_c(t, tau, sigma):
+    res = []
+
+    for t_i, tau_i in product(t, tau):
+        c = (ctypes.c_double * 3)(t_i, tau_i, sigma)
+        user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
+
+        area = quad(LowLevelCallable(lib.logGaussianCorrelation, user_data), -np.infty, np.infty)[0]
+
+        res.append(area)
+
+    res = np.asanyarray(res).reshape((t.shape[0], tau.shape[0]))
+
+    return res
+
+
 def _integrate_process_time(args):
     omega_i, tau_j, sigma = args
-    return quad(_integrand_time, -50, 50, args=(omega_i, tau_j, sigma))[0]
+    return quad(_integrand_time, -50, 50, args=(omega_i, tau_j, sigma), epsabs=1e-12, epsrel=1e-12)[0]
 
 
 def _integrand_time(u, t, tau, sigma):

src/nmreval/fit/minimizer.py

@@ -1,4 +1,3 @@
-import time
 import warnings
 from itertools import product
 

src/nmreval/io/read_old_nmr.py

@@ -9,11 +9,14 @@ import _compat_pickle
 
 from pickle import *
 from pickle import _Unframer, bytes_types, _Stop, _getattribute
+
+from nmreval.lib.logger import logger
+
 try:
     import bsddb3
     HAS_BSDDB3 = True
 except ImportError:
-    warnings.warn('bsdbb3 is not installed, reading legacy .nmr files is not possible.')
+    logger.warn('bsdbb3 is not installed, reading legacy .nmr files is not possible.')
     HAS_BSDDB3 = False
 
 

src/nmreval/models/fieldcycling.py

@@ -86,7 +86,7 @@ class EnergyFC(_AbstractFC):
     name = 'Energy distribution'
     params = ['C', 'T'] + EnergyBarriers.parameter
     bounds = [(0, None), (0, None), (0, None), (0, None)]
-    ralax = Relaxation(distribution=EnergyBarriers)
+    relax = Relaxation(distribution=EnergyBarriers)
 
 
 class _AbstractFCDipolar(_AbstractFC):

src/nmreval/nmr/coupling.py

@@ -14,7 +14,7 @@ from collections import OrderedDict
 from ..utils.constants import gamma_full, hbar_joule, pi, gamma, mu0
 
 
-__all__ = ['Quadrupolar', 'Czjzek', 'HeteroDipolar',
+__all__ = ['Coupling', 'Quadrupolar', 'Czjzek', 'HeteroDipolar',
            'HomoDipolar', 'Constant', 'CSA']