add loggausian to c func and move lookup of c library
This commit is contained in:
Dominik Demuth
parent
3e513a1231
commit
56b588293d
@ -1,15 +1,63 @@
|
||||
/* integrands used in quadrature integration with scipy's LowLevelCallables */
|
||||
#include <math.h>
|
||||
|
||||
/* 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);
|
||||
}
|
||||
|
||||
@ -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
|
||||
|
||||
27
src/nmreval/distributions/helper.py
Normal file
27
src/nmreval/distributions/helper.py
Normal file
@ -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')
|
||||
@ -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
|
||||
|
||||
@ -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()
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
|
||||
@ -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']
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user