From ead30d127a92968bb3efba51b29471526b19172e Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 10 Apr 2023 19:33:23 +0200 Subject: [PATCH 1/8] add C function for FFHS integration --- Makefile | 22 +++++++------- src/nmreval/clib/integrate.c | 15 ++++++++++ src/nmreval/distributions/intermolecular.py | 32 ++++++++++++++++++++- src/nmreval/fit/minimizer.py | 1 - 4 files changed, 57 insertions(+), 13 deletions(-) create mode 100644 src/nmreval/clib/integrate.c diff --git a/Makefile b/Makefile index 8c9b2d5..5f03e4f 100755 --- a/Makefile +++ b/Makefile @@ -2,37 +2,37 @@ #binaries PYUIC = pyuic5 -PYRCC = pyrcc5 +# PYRCC = pyrcc5 + +CC = gcc +CFLAGS = -O2 -fPIC +LDFLAGs = -shared + #Directory with ui files -RESOURCE_DIR = resources/_ui +RESOURCE_DIR = src/resources/_ui #Directory for compiled resources -COMPILED_DIR = nmreval/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) +ui : $(PYQT_UI) rcc: $(PNG_FILES) - $(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 clean: find . -name '*.pyc' -exec rm -f {} + diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c new file mode 100644 index 0000000..3cec66e --- /dev/null +++ b/src/nmreval/clib/integrate.c @@ -0,0 +1,15 @@ +/* integrands used in quadrature integration with scipy's LowLevelCallables */ + +/* integrand for distributions/intermolecular/FFHS */ +double ffhs_sd(double x, void *user_data) { + double *c = (double *)user_data; + double tau = c[1]; + double omega = c[0]; + + double u = x*x; + double res = u*u * c[1]; + + res /= 81. + 9.*u - 2.*u*u + u*u*u; + res /= u*u + omega*omega * tau*tau; + return res; +} diff --git a/src/nmreval/distributions/intermolecular.py b/src/nmreval/distributions/intermolecular.py index 1607361..c1b6354 100644 --- a/src/nmreval/distributions/intermolecular.py +++ b/src/nmreval/distributions/intermolecular.py @@ -1,8 +1,24 @@ +import ctypes +import os.path +from pathlib import Path + import numpy as np +from scipy import LowLevelCallable from scipy.integrate import quad 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') + class FFHS(Distribution): name = 'Intermolecular (FFHS)' @@ -24,7 +40,7 @@ class FFHS(Distribution): return ret_val @staticmethod - def specdens(omega, tau0, *args): + def specdens_py(omega, tau0, *args): 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 +49,17 @@ class FFHS(Distribution): return ret_val * 54 / np.pi + @staticmethod + def specdens_c(omega, tau0, *args): + 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) + 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 +75,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/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 From 3e513a12315c2dfed2f6773f930a159263327ae1 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Tue, 11 Apr 2023 19:06:35 +0200 Subject: [PATCH 2/8] update makefile --- Makefile | 17 +++++++++++------ src/nmreval/clib/integrate.so | Bin 0 -> 15624 bytes 2 files changed, 11 insertions(+), 6 deletions(-) create mode 100755 src/nmreval/clib/integrate.so diff --git a/Makefile b/Makefile index 5f03e4f..f0b77fa 100755 --- a/Makefile +++ b/Makefile @@ -4,11 +4,6 @@ PYUIC = pyuic5 # PYRCC = pyrcc5 -CC = gcc -CFLAGS = -O2 -fPIC -LDFLAGs = -shared - - #Directory with ui files RESOURCE_DIR = src/resources/_ui @@ -22,11 +17,21 @@ 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) +CC = gcc +CFLAGS = -O2 -fPIC +LDFLAGS = -shared + +C_DIR = src/nmreval/clib + all : ui ui : $(PYQT_UI) -rcc: $(PNG_FILES) +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 $@ diff --git a/src/nmreval/clib/integrate.so b/src/nmreval/clib/integrate.so new file mode 100755 index 0000000000000000000000000000000000000000..e726c62d1a9920fc871dd27baae10bc8c3c494f9 GIT binary patch literal 15624 zcmeHOTWl0n7(UyjP@u3CXeeN0B&c9yX;Ba*g)O}>U<+7E0tq-wySv>rd#SrqY4O1( z)dVCYeUa#!8ce_kB|fMxB-|tzH33aj5=|N-2^ude5k+HI|9|Fu+nMREF+PxJ&Pite z^IyLIob#QtJJaqt&$f4VR0RS8B`DSl)*^L|B_bUc*T@Qqh*&G8(|v(xuw$o8xF*#! z>4X$AFSk+6(Ng!zr(~BP;b*g*Qe!X`&Ayl9aj8POaZ|}R5_SSc)_dhEvc2O{0#`+0 zj2n}hR1$t4C??$k${@8YU{=&Dt#A8y`4 z%C)`4V)0snZWpRsNDu5qyO7^LL(X_OwD>vhl^DAo+I{+ZsILB(*+P)a*3d0P#t-Pa z6b~&c&^3NNw5vcm%lc34z}+8{&bIU9|B~cabdhb(XB2amu8*#Vp7?^UE95^%blN<7 zn6=wiH^qQrKrx^gPz)#r6a$I@#eiZ!F`yVw4BSQrmPG8vB0Z77qTsrvr0{Qpb%kP) z9U935$pXm{l3e(=UoL?I-6Akl7g$s~t!9jVfMo2QbY1+ZShRn-oL$#3yJ1Uc=Kh)y zv0>r5`yX7ogmH|gKkLOwio9<&RsG=B93@8Q^1BwB~fMP%~pcqgL zC>{b-Z^(sG$39yATa*j@H*OWPtnZVX zVzCdoL2icSeY9K-O8*Q;6?>(eDP@C{62_b2?J`Z;hol4T9+!67X>1qZe=j=yDdmy= z+a<7%-QOPe7jJ59eN~~8$z8{e9VX6e^zewm0x5uTG zu7H@#&y!^5b_>4zPuua6d47@X<@#S~T(8=BW{G1OReYsfp0_9_K#S8qU+>xWxutO- z)p=?6&&y@9SMkXT-?%=vfs@zUcXs^bbz89WPhJmGX~3GD%X_P@U|(6m{utTIowL|c zAzt+5th_ZCkB9rj#IXy*N*jHA@FGv=Goo3eKb7q@Q%2Ov=JJL)I3)VA>48+tibcbX zu6Cn+ zYa1UAY1_TErE7C5ndxW>deYQ|GpgEuF@; zj*cDeJw{K<#?E$1>peeWA3`aA(ByB9r{L6y??j6T=ZDjl*-O&OIWmD-CQD=N&kTlp z2a~DjieyySTEfgHL^wK}A*3U%oCxP)DU*#NY!5L!kg~%4Sxbbi*bv>Ba0<#sP0JMF zSVB^qh!PZ<19iA|5U5jLGo9?C_Oljere^Itz4^Qd(`=_{61^+)-)a0F!Ow2_*e^Y= z@Z946YrD%7w?cJ39`-9e%Yok~E0DnBIgj@Oq*+T9ezyP~?;nt3GLYM)k4))(2JaCO zvT<(U@tlOj`wNBGg>$%zP73D+9?w%qyubMSkM_F>%Q^A>RQf)W-%In!=PrK!-X^CT zuSp!FR7)wKL+!+^x-2iG332I!q6@{)bGRNy#NyLeTesuKP~a}4eS)n zEB)=E-;je8&RxX(RXR5)v)XfwI?GYI+<3f?Kn}Qrsg8DGK0+0KhGTu;{T$~C{qfKt zUnXqGg~$8NLH=M-L7bya=&-*|4qkg`AJ0{s+r;@y@D3b!<=>om{J`UVX`FVEa?$Mj z#0_}JH_7G3i-?SEbR<3}(pnl1@b8ktjmLAoLE;hT$3s6w6^=pwwTR= Date: Mon, 17 Apr 2023 17:31:34 +0200 Subject: [PATCH 3/8] 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'] From 4a9d50e8a4acf010fce82d429de75091e890a5f9 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 17 Apr 2023 17:43:18 +0200 Subject: [PATCH 4/8] refactor Log-Gaussian susceptibility --- src/nmreval/distributions/loggaussian.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index db50640..8c381a1 100644 --- a/src/nmreval/distributions/loggaussian.py +++ b/src/nmreval/distributions/loggaussian.py @@ -48,18 +48,17 @@ class LogGaussian(Distribution): _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 = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag) + res_imag = _integration_parallel(_omega, _tau, sigma, _integrate_process_real) + else: + res_real = None + res_imag = None - 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: ArrayLike, tau: ArrayLike, sigma: float): + def specdens(omega: ArrayLike, tau: ArrayLike, sigma: float) -> np.ndarray: _omega = np.atleast_1d(omega) _tau = np.atleast_1d(tau) @@ -106,8 +105,7 @@ def _integrate_specdens_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> n return res -def _integrate_process_imag(args): - omega_i, tau_j, sigma = args +def _integrate_process_imag(omega_i, tau_j, sigma): 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] From 2042148d0f80d466880223d96b66e2da48931db0 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Mon, 17 Apr 2023 20:17:54 +0200 Subject: [PATCH 5/8] C function for energy distribution spectral density --- src/nmreval/clib/integrate.c | 27 ++++++++++++++++- src/nmreval/distributions/energy.py | 45 ++++++++++++++++++++++------- 2 files changed, 60 insertions(+), 12 deletions(-) diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c index 4654530..1d221cf 100644 --- a/src/nmreval/clib/integrate.c +++ b/src/nmreval/clib/integrate.c @@ -1,6 +1,8 @@ /* integrands used in quadrature integration with scipy's LowLevelCallables */ #include +#define KB 8.617333262145179e-05 + /* FFHS functions */ double ffhsSD(double x, void *user_data) { double *c = (double *)user_data; @@ -18,7 +20,7 @@ double ffhsSD(double x, void *user_data) { /* 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; + return exp(- pow((log(tau/tau0) / sigma), 2) / 2.) / sqrt(2*M_PI)/sigma; } double logGaussianSD_high(double u, void *user_data) { @@ -61,3 +63,26 @@ double logGaussianCorrelation(double x, void *user_data) { return dist * exp(-t/uu); } + + +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); +} diff --git a/src/nmreval/distributions/energy.py b/src/nmreval/distributions/energy.py index a97a462..3ef11df 100644 --- a/src/nmreval/distributions/energy.py +++ b/src/nmreval/distributions/energy.py @@ -1,13 +1,18 @@ +import ctypes from itertools import product 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'] @@ -51,7 +56,7 @@ 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) + 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) - @@ -60,23 +65,41 @@ class EnergyBarriers(Distribution): return np.array(ret_val) @staticmethod - def specdens(omega, temperature, *args): + def specdens(omega, temperature, tau0: float, e_m: float, e_b: float): # 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) + if HAS_C_FUNCS: + ret_val = EnergyBarriers.spec_dens_c(omega, temperature, tau0, e_m, e_b) + else: + ret_val = EnergyBarriers.spec_dens_py(omega, 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 omega for tt in temperature]) + return ret_val.squeeze() - return ret_val + @staticmethod + def spec_dens_c(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 = (ctypes.c_double * 5)(o, tau0, e_m, e_b, t) + user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) + area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-10)[0] + res.append(area) + + return np.array(res) + + @staticmethod + def spec_dens_py(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray: + 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) + + e_axis = np.linspace(max(0., e_m-50*e_b), e_m+50*e_b, num=5001) + + ret_val = [simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis) for o in omega for tt in temperature] + + return np.array(ret_val) @staticmethod def mean(*args): From 2d472bd44e51a171420212896543980a7f8cee7f Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Wed, 19 Apr 2023 18:24:29 +0200 Subject: [PATCH 6/8] lower epsabs for integration of energy barrier spec dens --- Makefile | 4 ++-- src/nmreval/distributions/energy.py | 8 ++++---- src/nmreval/distributions/helper.py | 20 +++++++++++--------- src/nmreval/models/fieldcycling.py | 2 +- 4 files changed, 18 insertions(+), 16 deletions(-) diff --git a/Makefile b/Makefile index f0b77fa..24001f8 100755 --- a/Makefile +++ b/Makefile @@ -23,8 +23,8 @@ LDFLAGS = -shared C_DIR = src/nmreval/clib -all : ui - +all : ui compile + ui : $(PYQT_UI) rcc : $(PNG_FILES) diff --git a/src/nmreval/distributions/energy.py b/src/nmreval/distributions/energy.py index 3ef11df..1120b78 100644 --- a/src/nmreval/distributions/energy.py +++ b/src/nmreval/distributions/energy.py @@ -1,5 +1,5 @@ -import ctypes from itertools import product +from ctypes import c_double, cast, pointer, c_void_p import numpy as np from scipy import LowLevelCallable @@ -82,9 +82,9 @@ class EnergyBarriers(Distribution): def spec_dens_c(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 = (ctypes.c_double * 5)(o, tau0, e_m, e_b, t) - user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p) - area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-10)[0] + c = (c_double * 5)(o, tau0, e_m, e_b, t) + user_data = cast(pointer(c), c_void_p) + area = quad(LowLevelCallable(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-13)[0] res.append(area) return np.array(res) diff --git a/src/nmreval/distributions/helper.py b/src/nmreval/distributions/helper.py index 27c4db2..11c0ba9 100644 --- a/src/nmreval/distributions/helper.py +++ b/src/nmreval/distributions/helper.py @@ -1,24 +1,26 @@ from pathlib import Path -import ctypes +from ctypes import CDLL, c_double, c_void_p from ..lib.logger import logger - lib = None try: - lib = ctypes.CDLL(str(Path(__file__).parents[1] / 'clib' / 'integrate.so')) + lib = 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) + lib.ffhsSD.restype = c_double + lib.ffhsSD.argtypes = (c_double, 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) + lib.logGaussianSD_high.restype = c_double + lib.logGaussianSD_high.argtypes = (c_double, c_void_p) + lib.logGaussianSD_low.restype = c_double + lib.logGaussianSD_low.argtypes = (c_double, c_void_p) + + lib.energyDist_SD.restype = c_double + lib.energyDist_SD.argtypes = (c_double, c_void_p) HAS_C_FUNCS = True logger.info('Use C functions') 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): From dd1c26e285f602faf14d28d640ae1502add831cd Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 23 Apr 2023 19:55:53 +0200 Subject: [PATCH 7/8] more work on loggaussian --- AppImageBuilder.yml | 23 +++++++++++------ src/nmreval/clib/integrate.c | 31 +++++++++++++++++++++-- src/nmreval/clib/integrate.so | Bin 15624 -> 16088 bytes src/nmreval/distributions/loggaussian.py | 22 ++++++---------- 4 files changed, 52 insertions(+), 24 deletions(-) diff --git a/AppImageBuilder.yml b/AppImageBuilder.yml index dae00b6..ec7aec7 100644 --- a/AppImageBuilder.yml +++ b/AppImageBuilder.yml @@ -34,10 +34,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 @@ -55,7 +55,6 @@ AppDir: - libqt5test5 - libqt5xml5 - qtbase5-dev-tools - - qtchooser - pyqt5-dev-tools - qtchooser - libavahi-client3 @@ -63,10 +62,18 @@ AppDir: - 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/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c index 1d221cf..f23b5b2 100644 --- a/src/nmreval/clib/integrate.c +++ b/src/nmreval/clib/integrate.c @@ -23,7 +23,7 @@ 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 logGaussian_imag_high(double u, void *user_data) { double *c = (double *)user_data; double omega = c[0]; @@ -36,7 +36,7 @@ double logGaussianSD_high(double u, void *user_data) { return dist * omega * uu / (pow(uu, 2) + pow(omega, 2)); } -double logGaussianSD_low(double u, void *user_data) { +double logGaussian_imag_low(double u, void *user_data) { double *c = (double *)user_data; double omega = c[0]; @@ -50,6 +50,33 @@ double logGaussianSD_low(double u, void *user_data) { 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; diff --git a/src/nmreval/clib/integrate.so b/src/nmreval/clib/integrate.so index e726c62d1a9920fc871dd27baae10bc8c3c494f9..2580e2490d0ad56587711645b604a1d52cf42f0b 100755 GIT binary patch literal 16088 zcmeHOdu$ZP8K3im953IILSi1Rx5Nm8$YER&Xd__3XLD(IIADs6fc4pTekA8Rdbb9f zqT+-m3PMPh`bed!DpC7KR7HhWM4_rw8H_@j(iST9omN*$q;o}EhqNiQf$Q&^neV)} zoRQi;s!H9FX21D<^S$Pq-SwPz=4jKF<}#m8a48q-1#u%mg_yddZLMyAm}0G%j_U$3 zU&>ZYcHP^WrXi-(F{khqWG4zgDLhsjk)b+_YB9F(p|lE$`SdYMTae_tL!z7ED?1 zDd=&Ewa9%5r+tU^H&tHF=Pr->;p%bbS2fVE*0^eDdIu$U33n6ZSg>8P_&s|TN4c6-I|)}C4a)gu=V-sul-7(e@m_txF% zJooUPXKp-K2azeu)9B%`IE9+Q7*0U9S4@GE4MWQ~NfJyc`!pWIDX81u2Ypb?5C$a& z<%A+44N7?MK?TUmXFsIy21VrjlKki#1vF{=B~WQc%adIyyHE*l(em9A72;bO zAJOugH2yr}I-UWI3(M;1PbIC49ZlPoC9Ee}w_34yI^NTlvE%91t&NFPGTs{POvF`R z=^m?VFlu%8C8LSH!*S8w-J97N7Kv2Pwp6-5nh3)daUwo=favDvKqk`{P3{a^y?s5s zE+~;Y2y31wfj;@k(@M_j@oe%PsYI*)_kiyn0WVABG^D+W|#!WG===L zQhxPeBFpYZj!#6+|9L~?{H^kc?~TaSyY^fdxPzbXWWb-#b_Z(k(}4Ad*KLHJ7+AS0 za%|m?=|_I#_&fHj$cc4}f#>!i#9S{fZ&c0)@7u?EOvgWi&hAj4=5rL47`Yi}AH5j} z&ih7{5TGePE&+f)K)uo(sGC5|-VC%&fG0Zdm=tXP5qQEbqW^8s-=YHB`*X-#K)rl3 zu>Y5+ze4}m*rY>1=P`%XHv@B8YyG)>ck_8b;ad1@AMFn8u6$D2uY3eF_gU3l`4I_{ z7-qYAtb5WPY}Fy;rKsB}_T51;GEK_QiPe&yEww@IZKX69gkwkoMr?>NSa}J*NFcJP z7NLygZr`Dpf>JrSs1_dKq5UH2b<{D`-(Y|k^!~B25$(SQ3feym?REKmU0E#ng$Ra{ z8&R;7hC4sP`5~+$q;{mQozjYFXq#P!@Kh9>&U?mxfB$QgBbf*_diW3Kg2nzRb05Au zOQu3j$o{u!@^ZI7pnSBew9ZU;mY zN#ZVhbZ;Pz*y#jh2wB=?#{{iu!bMc>pmiO}c6`-d;!=PEvz4RXi6#W@h_zv$7zQLex;j9Pz zBj~F|-y+lw`0M}`Rl(+}>Q4q{9`X;14GSJ${;}G{1e1Ok^~;dsZ+e8&H(iRm1+ntC zQj$G@`YJH_y){y09I5h48!rOC4*V#n2v^mYbyoRF{uJsKj7xq$CA?SRH14y&<@Z^_ z`xSl#_#%vZ6R3dw0c}4AwiS2-;o&L+;%e*|>k;7cJ2KgM0UciPGT>#v%Yc^wF9Ti% zybO35@G{_Kz{|k@Kn8e!CGV%iM> zoQP_dY#^ZZ5NrQ)AzTrW#IRpKlK>?)<~oQ8}KgL2Oo{ol!l z{JZmfQRD7;xeC0D-e#`6{Q@|>txZEV3YR~kRfe~w`@Wq3?i2Uk57RJUmuYEh25>sB z;6C>d(75}&!UEui)GWR%)AFNw|Gxac9Yoe>+|Vy+6zA^>=I`C3?_#W_5=Qb`(VMlbO?TkHu9XnlO@>qst_gh_bFe8KW8?lttlSp+& z6IRSlr88D^U{G|W`VS=Hc03lUFD*fbHTtY*IvpLd;z>I_B)Zel{9si>V6AvzR8pl%&V*dioz8$vVGlLBkU2XUpDhEOUNwWA^w z@6~hL8-o$%l&YF!rNM1@i1zn&A(E6$mJxuo)0xSL5N@-6+&#sQ=zqk~I*`^#iLqSkMg?C z1V*-5pVuQydHo3yS;&U7M8J@3*3aqzGUfH8JAU@RL+fwSMtOb8bWr1@Pp>>4KkL)C z5sHb|&xMzleu?qA7b3Fg@R$|5(Q@h=T90W*K~gpN-N>}lr9Z4Ym~y;qhIyvXxb%6w zz;vN2VD`)QKJC)yzjsWXXR%1fWlkK@`aFNaASo0>UWKVh;p+BJ-0;$_y_Ej5Ym~yzoJP@h>;DU4zbDK9 delta 1980 zcmZ`)Z)jUp6u&pGX`A%lOVjZ3+StLHLUEGYVsI;up zI+Q&eLgp-{9+-bV>8P?QNdJIRL9qKHg?-TJ2iGrBKM4J12a4Me;<@kLOP7Eb?)}|+ z&hMOe?z!jQcjwLgxq`PdEfce)Eg+ju@(Dtk#E1h&m3S9Yxzd%4O{*x9&Uf!h_hoxl;!=Lo}eP?t5!cR99O zV$u%lbXfntNDS`LfG(Rav87!5=h)>TXS#U$?{{7s(~f+7^E>8-#bw9N$ZTIMDF${X zW9|CVzl(62w!$4-2ByOmkdT^OM^Vv;xg;)=LvrA#IGLF^bDAwhMyf3x+2xtb^BoKM zGM9VP7E#E(pxfQ5NFUSfq$~!D&jv$w1)j0jI_?SSLcg~5lm27P`*LXTivS_Hf<-4J zGoRc-L1y0mGlSlHv!2W)`(-wN`}LXCs3K7|aZETvk_{VH*=2G%0T+E0Ihw_H7x ztbN(^Y|WZnIUqg0FO8;HnWJgiA{}hUsf1nLEf&*=8sa44E}}u;W39u!`Dd{x!%gkc zN0xCWQidF}&c+ZYk(^uG(BG@YPmEHhNdK}7C;V-$Nq%H1v%zM-MZb^cVaXrRz3k)I zSe{wSv*N(!w1AAmeSfF3im&D*I1%`q7Jvp7>V^K`d33W8R06^c=b0d3LCtaZrSeMf zbt_U8(LnE3Xb?w_gtrt9cPpw0w$!4A5LiK0Es;&EBQZpU8t^+x8%Q2e!MJdxy2v&! zlTsgzkP0c1n_%1_g=#o@Y#st9!~0=5RM~<7@ahz1^*FI%$xp-WkO%dmdWAlU6H1zJ z8#L?*V-URBh5kiuFOEaVUp4LO@f-%jhom~Ph+A>PVz^n7NdY#(o)$OJ(t>T_4bz4N ze;4(5T4k~RYSLJG^d67uB;{$`HSHT=xW@ydktXVa!x4oB*r+~@7q8@MTc0FHn09SB zQc3G4CQi*|CvsD-WM7&(_dEon(W+c-@(f;on2DZ(-=g2b2eB+|hx@TWqo^T)1p@J2 zzf07Jz;Vci-SBR2t5Z}7me*?UW6&B;&{a4d_rVRD3|;YtW^&97f3Hk~eEUIIiFfYH zmTAR73hbcpAl|0GP(~9u66k>Sn9uPtXEj3hpkhX_8JC;Saw)$aWbb{3)mTT@g)((a zLjg8{)7BEJnD#2_Jg!)AiOV6LPXoFWo$5!ZX?RF0$|}#a`e&$NmsQw`pNaU%jKF%r zN5^0*k)mILn)Km4U61?JZ_tK4u)qaUpxq{FIhx!gUKJ2k1&%x_`QOARbwYtU6y!Jmy`Zxyw diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index 8c381a1..7f8987b 100644 --- a/src/nmreval/distributions/loggaussian.py +++ b/src/nmreval/distributions/loggaussian.py @@ -52,8 +52,8 @@ class LogGaussian(Distribution): res_real = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag) res_imag = _integration_parallel(_omega, _tau, sigma, _integrate_process_real) else: - res_real = None - res_imag = None + res_real = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag) + res_imag = _integration_parallel(_omega, _tau, sigma, _integrate_process_real) return (res_real + 1j * res_imag).squeeze() @@ -95,8 +95,8 @@ def _integrate_specdens_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> n 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] + area = quad(LowLevelCallable(lib.logGaussianSD_high, user_data), 0, np.infty, epsabs=1e-12, epsrel=1e-12)[0] + area += quad(LowLevelCallable(lib.logGaussianSD_low, user_data), -np.infty, 0, epsabs=1e-12, epsrel=1e-12)[0] res.append(area) @@ -105,9 +105,10 @@ def _integrate_specdens_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> n return res -def _integrate_process_imag(omega_i, tau_j, sigma): - 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] +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), 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 @@ -166,10 +167,3 @@ 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() From 8086dd5276be4ec52115b72900ef5b63bf94f5d1 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sat, 29 Apr 2023 20:40:32 +0200 Subject: [PATCH 8/8] finalized c functions --- src/nmreval/clib/integrate.c | 47 +++++++- src/nmreval/clib/integrate.so | Bin 16088 -> 16344 bytes src/nmreval/distributions/energy.py | 133 +++++++++++------------ src/nmreval/distributions/helper.py | 29 ++++- src/nmreval/distributions/loggaussian.py | 22 ++-- 5 files changed, 146 insertions(+), 85 deletions(-) diff --git a/src/nmreval/clib/integrate.c b/src/nmreval/clib/integrate.c index f23b5b2..28c591a 100644 --- a/src/nmreval/clib/integrate.c +++ b/src/nmreval/clib/integrate.c @@ -1,7 +1,7 @@ /* integrands used in quadrature integration with scipy's LowLevelCallables */ #include -#define KB 8.617333262145179e-05 +const double KB = 8.617333262145179e-05; /* FFHS functions */ double ffhsSD(double x, void *user_data) { @@ -91,7 +91,7 @@ double logGaussianCorrelation(double x, void *user_data) { 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; } @@ -113,3 +113,46 @@ double energyDist_SD(double x, void *user_data) { 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 index 2580e2490d0ad56587711645b604a1d52cf42f0b..9a8f34467ada92cbe20305d4b174125247930bce 100755 GIT binary patch delta 2880 zcma)8eN0nV6o0Q!s0E=dNFmWtv*ts?XcjjY z1+Q6_FsJb^On0Rc^5yr7-%1_x;iEl7;Ui%|JhBL-}uH)BZ8D%!{G!(Z8 zTG_&msyhO&S;w2}zRq2`^+w)%r30RkNWPS6YnY&v{)?D@Ut{G@3DgNv5A9rUAT$tR}er)5(D<2%OF1g-$ zd~frp0|)-h+}Mnz?3|j$Lb~>VOH~xJ60ayWu=p7SD7SzXiZntH+azy-J@L*ILF|zH z4IuHMWI^OgevGsaN&7D(Um@(-o+$IMK5?(yW=Eqk%|ZFJDM^?%FmF`4=c;6zrTZ`ii{3K3s=mQnTUA$ITc`M%YUg-YuUb{_Z4|7p zY3&PSOO4m}vMI4_Uv$yTn#P)z+Vv&%tNib-UbXxk(8`O{T!ea#dZVeOrN-y=*EcmP z^JcM2(Lb}x`gnFkAD?EHnhBh4@j&|@AN)-pHLT0DSXB=)k2y=-#@ftPDeosJN`In9 zKM@$-_*B#S6D#%O*gxiMHHBqb9O@ai$nuryWcMu-Vm$h|;V6UMmefG@0GXPj?cuTf z&a}{O=ub4wx0}Mzf~c^kg0S05H}_^C!X=E}{=~(X5um3hG2N|{ZoY{L(_6Hz+=NRe z9L}dz=4x_gUI=F=o=_ZkCj&8A0y=N6;4AE=6^U7Wm@A8 zR7h6(wQ=l3QpPC%^Gi=+V)zm;BqaRFW2R3oR~5E2*&R6Zh!%UId6uBDgjC<03d1mT zAR(r;q%%OO%59)eLy(VapF`5C?OC|`?KWIn?Wy4=?$3lThmKm^-mbTU>pxG z;)VfJ1>ZKf0o#W)`V{7)>N_?K6R=xGHEfae*W4 zc=qv4^&`z^AM#`o3wfphfKMXjr0++J3y1_IJUVS9wv=xnipLyphr#J^C-)M`TD#2H z{}WWfGv?zN^H(Be#X^#wr;PI?2EkAkB!?nm^oGX?LpY$~Fbo?k4iq1f;OPmE5yFsc zN?}j%VDf&0Vb7W0cRl?)rFf?;D{@eVSe9*MfDU&$lF#1Kw7o#*ZB3(NO}V3K*8t7H zd%zGd9!H|PU(@Ws5O6xsc~{dG0`~%aK-)b{+X8F`ehEATJPdT-*R)f>YryL?|F@<+ zK>>E!#;RGY+m@NO3<@IJ5Ku8;&|#a()QqGwKXkoGbZLUEHNl|JX-Y@sK{h$V8c2pZ zeo&SqICXUi2AVyB(F&d)ti+cK{w(-&;2VGn@yi4s#^WTsRZE0#5WE7lyTJ2<8_`$G z=&9gOfk(rrCr38eb6HMI+{47Sh@+2%qf+obq^VQzZe5$;mytbq-jPVqr-ENc_DCP^ zJS7Qk9ZN6-yEAUau`)-NZn2eZa8!;CMUyv0i4mam$(By7mVVP!_Q(07R9*FG zwsXpM@ELhlGAhcOZGP4SrK@yyC(jzl$1A<6*%pC+RWqLG{)aSe zKX@CuMtQzj-ORf4taCf0<J|1jt`F#{OV+XU zoDum-Dt3Z2wnrMa%)~^=m&q=E$nJRie=EAz6hxOv)IBhdj566 zCKVPBF`!R|;L5`7s#?t+PR%w)N-=^Dn9G%Ii4<%EQJ@=K4&$=|j$vJ{0%N3J_s0bL zT`{VzFNQrTau^p|MY2$2xzlE&;?_^gHox?}%8r7)XcMLpY!5!1R;*LcGS|#Pqq{&% XzOG{%W)>Qzil_6;J0*#I=8FFpe`F69 delta 1481 zcmZ8heN0HXsOIr(@m01*JbE{G|Vo_Amd}gnA!+`ePS%ldy1ogfB7FLYhnr&=dCaltPrckcQM`GQ{_(Oa=2X0gF{#i@ z5-yQ-8`=#NKl&;ELM`^Ld}q4V(Sm@TZ&KAftw@Z}3&P6SBj3{bR7MUnV8+y+J6WFRcwfEEojQSp2t9^cOqH4Bt-3XV78s zkZ~BZ9N4nU#4vethrR|e_u9Xbls9)_v3nLDF+i==ORm8S){lt?{<7A?dFz&FXoG6_ z-EJP|Hn%Zvw#40hj#ucU+m7^0Zob>RzjCRR*2XZ3x_;_5dlqvAXtOm{j>&-zeJ5gW zZe1=Ym+I2*DTU8%uFCj>so!CGSv@jXVjLc3lY9#JTQUhWJF^-r;Ou~w@mb?fBoIe1UODIoZ)o#q9gWiB6Xw&elOK{x6hheWsH91F2M!9AV`8MQ(ir$HGnWq={ z6!K>k-miPpWW+uAc@y}Hzk!6;qnozDNpDv*ueK7s>~u`&6r|fj#}Wf;-sh|ndc~!E zKu4ApJAK`ROu>R*2x>)?m@*eJdQ?5gXuN@$&;u7ZMqgnp;3YBmG9Zu)&gJ5bzzosN zy5PO7=a6?b3RDzn+*xeT2Jiz45lW`$J#en+bO4oEi& zfp#^@t6e7_0xAZzwXBMKM9q4gu0=kS=dU2uC^)`Wl?ioS?fAc?8Zp+y5%W;$=*!H& z;81_6FWx(R6vR-2_5V#w`x?8rMz=2xj&M{w9Daw8I24+F7VUI2?1AQp&!TO=hD(sz zx*6V%RB+k>WLcPrG;!J~PUyva#6Wb9vat4qmz!4qPq@{=mYuloCt7@#m{a*m4e`H; kY(`}@w2Ot7s18$_qap6%HYH_+z{zOH7*f0+)VC@951S8e7ytkO diff --git a/src/nmreval/distributions/energy.py b/src/nmreval/distributions/energy.py index 1120b78..40f4c7f 100644 --- a/src/nmreval/distributions/energy.py +++ b/src/nmreval/distributions/energy.py @@ -28,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 @@ -56,88 +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, tau0: float, e_m: float, e_b: float): - # in contrast to other spectral densities, it's omega and temperature - - omega = np.atleast_1d(omega) - temperature = np.atleast_1d(temperature) - if HAS_C_FUNCS: - ret_val = EnergyBarriers.spec_dens_c(omega, temperature, tau0, e_m, e_b) + 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 = EnergyBarriers.spec_dens_py(omega, temperature, tau0, e_m, e_b) + 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.squeeze() + return ret_val @staticmethod - def spec_dens_c(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(lib.energyDist_SD, user_data), 0, np.infty, epsabs=1e-13)[0] - res.append(area) - - return np.array(res) + def mean(temperature, tau0, ea): + return tau0*np.exp(ea/(kB*temperature)) @staticmethod - def spec_dens_py(omega: np.ndarray, temperature: np.ndarray, tau0: float, e_m: float, e_b: float) -> np.ndarray: - 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) - - e_axis = np.linspace(max(0., e_m-50*e_b), e_m+50*e_b, num=5001) - - ret_val = [simpson(integrand(e_axis, o, tau0, e_m, e_b, tt), e_axis) for o in omega for tt in temperature] - - return np.array(ret_val) - - @staticmethod - def mean(*args): - return args[1]*np.exp(args[2]/(kB*args[0])) - - @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 index 11c0ba9..b165f27 100644 --- a/src/nmreval/distributions/helper.py +++ b/src/nmreval/distributions/helper.py @@ -14,16 +14,35 @@ try: lib.ffhsSD.argtypes = (c_double, c_void_p) # Log-Gaussian integrands - lib.logGaussianSD_high.restype = c_double - lib.logGaussianSD_high.argtypes = (c_double, c_void_p) - lib.logGaussianSD_low.restype = c_double - lib.logGaussianSD_low.argtypes = (c_double, c_void_p) + 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') \ No newline at end of file + logger.info('Use python functions') + diff --git a/src/nmreval/distributions/loggaussian.py b/src/nmreval/distributions/loggaussian.py index 7f8987b..10c3615 100644 --- a/src/nmreval/distributions/loggaussian.py +++ b/src/nmreval/distributions/loggaussian.py @@ -49,8 +49,8 @@ class LogGaussian(Distribution): _tau = np.atleast_1d(tau0) if HAS_C_FUNCS: - res_real = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag) - res_imag = _integration_parallel(_omega, _tau, sigma, _integrate_process_real) + 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) @@ -63,7 +63,7 @@ class LogGaussian(Distribution): _tau = np.atleast_1d(tau) if HAS_C_FUNCS: - ret_val = _integrate_specdens_c(_omega, _tau, sigma) + ret_val = _integrate_susc_imag_c(_omega, _tau, sigma) else: ret_val = _integration_parallel(_omega, _tau, sigma, _integrate_process_imag) @@ -88,15 +88,23 @@ def _integration_parallel(x: np.ndarray, tau: np.ndarray, sigma: float, func: Ca return res -def _integrate_specdens_c(omega: np.ndarray, tau: np.ndarray, sigma: float) -> np.ndarray: +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(lib.logGaussianSD_high, user_data), 0, np.infty, epsabs=1e-12, epsrel=1e-12)[0] - area += quad(LowLevelCallable(lib.logGaussianSD_low, user_data), -np.infty, 0, epsabs=1e-12, epsrel=1e-12)[0] + 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) @@ -139,7 +147,7 @@ def _integrate_correlation_c(t, tau, sigma): 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):