From 56b588293d9536ba08d00db69f07a3a81abf4598 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 17 Apr 2023 17:31:34 +0200 Subject: [PATCH] add loggausian to c func and move lookup of c library --- src/nmreval/clib/integrate.c | 56 +++++++++++- src/nmreval/distributions/base.py | 2 +- src/nmreval/distributions/helper.py | 27 ++++++ src/nmreval/distributions/intermolecular.py | 21 ++--- src/nmreval/distributions/loggaussian.py | 94 ++++++++++++++++----- src/nmreval/io/read_old_nmr.py | 5 +- src/nmreval/nmr/coupling.py | 2 +- 7 files changed, 166 insertions(+), 41 deletions(-) create mode 100644 src/nmreval/distributions/helper.py diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c index 3cec66e..4654530 100644 --- a/src/nmreval/clib/integrate.c +++ b/src/nmreval/clib/integrate.c @@ -1,15 +1,63 @@ /* integrands used in quadrature integration with scipy's LowLevelCallables */ +#include -/* integrand for distributions/intermolecular/FFHS */ -double ffhs_sd(double x, void *user_data) { +/* FFHS functions */ +double ffhsSD(double x, void *user_data) { double *c = (double *)user_data; - double tau = c[1]; + double omega = c[0]; + double tau = c[1]; double u = x*x; - double res = u*u * c[1]; + 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 logGaussianSD_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 logGaussianSD_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 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); +} 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/helper.py b/src/nmreval/distributions/helper.py new file mode 100644 index 0000000..27c4db2 --- /dev/null +++ b/src/nmreval/distributions/helper.py @@ -0,0 +1,27 @@ + +from pathlib import Path +import ctypes + +from ..lib.logger import logger + + + +lib = None +try: + lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so')) + + # FFHS integrand + lib.ffhsSD.restype = ctypes.c_double + lib.ffhsSD.argtypes = (ctypes.c_double, ctypes.c_void_p) + + # Log-Gaussian integrands + lib.logGaussianSD_high.restype = ctypes.c_double + lib.logGaussianSD_high.argtypes = (ctypes.c_double, ctypes.c_void_p) + lib.logGaussianSD_low.restype = ctypes.c_double + lib.logGaussianSD_low.argtypes = (ctypes.c_double, ctypes.c_void_p) + + HAS_C_FUNCS = True + logger.info('Use C functions') +except OSError: + HAS_C_FUNCS = False + logger.info('Use python functions') \ No newline at end of file diff --git a/src/nmreval/distributions/intermolecular.py b/src/nmreval/distributions/intermolecular.py index c1b6354..6ed193e 100644 --- a/src/nmreval/distributions/intermolecular.py +++ b/src/nmreval/distributions/intermolecular.py @@ -1,24 +1,15 @@ import ctypes -import os.path -from pathlib import Path import numpy as np from scipy import LowLevelCallable from scipy.integrate import quad +from .helper import HAS_C_FUNCS, lib from .base import Distribution -try: - print(str(Path(__file__).parents[1] / 'clib' / 'integrate.so')) - lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so')) - lib.ffhs_sd.restype = ctypes.c_double - lib.ffhs_sd.argtypes = (ctypes.c_double, ctypes.c_void_p) - HAS_C_FUNCS = True - print('Use C functions') -except OSError: - HAS_C_FUNCS = False - print('Use python functions') +# 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)' @@ -40,7 +31,7 @@ class FFHS(Distribution): return ret_val @staticmethod - def specdens_py(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) @@ -50,12 +41,12 @@ class FFHS(Distribution): return ret_val * 54 / np.pi @staticmethod - def specdens_c(omega, tau0, *args): + 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.ffhs_sd, user_data) + func = LowLevelCallable(lib.ffhsSD, user_data) res.append(quad(func, 0, np.infty)[0]) return np.array(res) * 54 / np.pi diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index e1736c8..db50640 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,14 +36,12 @@ 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): @@ -54,23 +59,51 @@ class LogGaussian(Distribution): return ret_val.squeeze() @staticmethod - def specdens(omega, tau0, sigma): + def specdens(omega: ArrayLike, tau: ArrayLike, sigma: float): _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_specdens_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_specdens_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> np.ndarray: + 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(lib.logGaussianSD_high, user_data), 0, np.infty)[0] + area += quad(LowLevelCallable(lib.logGaussianSD_low, user_data), -np.infty, 0)[0] + + res.append(area) + + res = np.asanyarray(res).reshape((omega.shape[0], tau.shape[0])) + + return res def _integrate_process_imag(args): @@ -89,6 +122,22 @@ 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] @@ -119,3 +168,10 @@ def _integrand_freq_real_low(u, omega, tau, sigma): def _integrand_freq_real_high(u, omega, tau, sigma): uu = np.exp(-2*u) return LogGaussian.distribution(np.exp(u), tau, sigma) * uu / (uu + omega**2) + + +if __name__ == '__main__': + import matplotlib.pyplot as plt + xx = np.logspace(-5, 5) + plt.loglog(xx, LogGaussian.specdens(xx, 1, 2)) + plt.show() 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/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']