From 8086dd5276be4ec52115b72900ef5b63bf94f5d1 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sat, 29 Apr 2023 20:40:32 +0200 Subject: [PATCH] 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):