IPKM/nmreval
17
2
Fork 1
Code Issues 42 Pull Requests Actions Packages Projects Releases Wiki Activity

more work on loggaussian

Browse Source
This commit is contained in:
Dominik Demuth
2023-04-23 19:55:53 +02:00
parent 2d472bd44e
commit dd1c26e285
4 changed files with 52 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
src/nmreval/clib/integrate.c
View File

@ -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;

BIN
src/nmreval/clib/integrate.so
View File

Binary file not shown.

22
src/nmreval/distributions/loggaussian.py
View File

@ -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()