diff --git a/AppImageBuilder.yml b/AppImageBuilder.yml index 8bcd034..80b4572 100644 --- a/AppImageBuilder.yml +++ b/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}' diff --git a/Makefile b/Makefile index 78133ec..2fc5214 100755 --- a/Makefile +++ b/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 diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c new file mode 100644 index 0000000..28c591a --- /dev/null +++ b/src/nmreval/clib/integrate.c @@ -0,0 +1,158 @@ +/* integrands used in quadrature integration with scipy's LowLevelCallables */ +#include + +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); +} + diff --git a/src/nmreval/clib/integrate.so b/src/nmreval/clib/integrate.so new file mode 100755 index 0000000..9a8f344 Binary files /dev/null and b/src/nmreval/clib/integrate.so differ diff --git a/src/nmreval/distributions/base.py b/src/nmreval/distributions/base.py index ae9b72e..68b10bf 100644 --- a/src/nmreval/distributions/base.py +++ b/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 diff --git a/src/nmreval/distributions/energy.py b/src/nmreval/distributions/energy.py index a97a462..40f4c7f 100644 --- a/src/nmreval/distributions/energy.py +++ b/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) diff --git a/src/nmreval/distributions/helper.py b/src/nmreval/distributions/helper.py new file mode 100644 index 0000000..b165f27 --- /dev/null +++ b/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') + diff --git a/src/nmreval/distributions/intermolecular.py b/src/nmreval/distributions/intermolecular.py index 1607361..6ed193e 100644 --- a/src/nmreval/distributions/intermolecular.py +++ b/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 diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index e1736c8..10c3615 100644 --- a/src/nmreval/distributions/loggaussian.py +++ b/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): diff --git a/src/nmreval/fit/minimizer.py b/src/nmreval/fit/minimizer.py index 2655e4a..57e1b56 100644 --- a/src/nmreval/fit/minimizer.py +++ b/src/nmreval/fit/minimizer.py @@ -1,4 +1,3 @@ -import time import warnings from itertools import product diff --git a/src/nmreval/io/read_old_nmr.py b/src/nmreval/io/read_old_nmr.py index fd8f1e4..83c283a 100644 --- a/src/nmreval/io/read_old_nmr.py +++ b/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 diff --git a/src/nmreval/models/fieldcycling.py b/src/nmreval/models/fieldcycling.py index 1da10c5..4d5221e 100644 --- a/src/nmreval/models/fieldcycling.py +++ b/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): diff --git a/src/nmreval/nmr/coupling.py b/src/nmreval/nmr/coupling.py index 5e45e54..b3534dc 100644 --- a/src/nmreval/nmr/coupling.py +++ b/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']