From dd1c26e285f602faf14d28d640ae1502add831cd Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 23 Apr 2023 19:55:53 +0200 Subject: [PATCH] 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()