finalized c functions
This commit is contained in:
Dominik Demuth
parent
dd1c26e285
commit
8086dd5276
@ -1,7 +1,7 @@
|
||||
/* integrands used in quadrature integration with scipy's LowLevelCallables */
|
||||
#include <math.h>
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
|
||||
Binary file not shown.
@ -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)
|
||||
|
||||
@ -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')
|
||||
logger.info('Use python functions')
|
||||
|
||||
|
||||
@ -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):
|
||||
|
||||
Loading…
Reference in New Issue
Block a user