From dab5eac5490e0b039348a878ea487e190ee96368 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Mon, 24 Aug 2020 16:36:39 +0100 Subject: [PATCH 01/25] mbar overlap matrix --- docs/visualisation.rst | 28 +++++++ ...hemlyb.visualisation.plot_mbar_omatrix.rst | 16 ++++ src/alchemlyb/estimators/mbar_.py | 4 + src/alchemlyb/tests/test_visualisation.py | 17 ++++ src/alchemlyb/visualisation/__init__.py | 1 + src/alchemlyb/visualisation/mbar_martix.py | 78 +++++++++++++++++++ 6 files changed, 144 insertions(+) create mode 100644 docs/visualisation.rst create mode 100644 docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst create mode 100644 src/alchemlyb/tests/test_visualisation.py create mode 100644 src/alchemlyb/visualisation/__init__.py create mode 100644 src/alchemlyb/visualisation/mbar_martix.py diff --git a/docs/visualisation.rst b/docs/visualisation.rst new file mode 100644 index 00000000..c84cfcae --- /dev/null +++ b/docs/visualisation.rst @@ -0,0 +1,28 @@ +Visualisation of the results +============================ +It is quite often that the user want to visualise the results to gain +confidence on the computed free energy. **alchemlyb** provides various +visualisation tools to help user to judge the estimate. + +.. _plot_overlap_matrix: +Overlap Matrix of the MBAR +-------------------------- +The accuracy of the :class:`~alchemlyb.estimators.MBAR` estimator depends on +the overlap between different lambda states. The overlap matrix from the +:class:`~alchemlyb.estimators.MBAR` estimator could be plotted to check +the degree of overlap. It is recommended that there should be at least +**0.05** overlap between neighboring states. :: + + >>> import pandas as pd + >>> from alchemtest.gmx import load_benzene + >>> from alchemlyb.parsing.gmx import extract_u_nk + >>> from alchemlyb.estimators import MBAR + + >>> bz = load_benzene().data + >>> u_nk_coul = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']]) + >>> mbar_coul = MBAR() + >>> mbar_coul.fit(u_nk_coul) + + >>> from alchemlyb.visualisation.mbar_martix import plot_mbar_omatrix + >>> ax = plot_mbar_omatrix(mbar_coul.overlap_maxtrix) + >>> ax.figure.savefig('O_MBAR.pdf', bbox_inches='tight', pad_inches=0.0) diff --git a/docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst b/docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst new file mode 100644 index 00000000..ccdaa6af --- /dev/null +++ b/docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst @@ -0,0 +1,16 @@ +Plot Overlap Matrix from MBar +============================= + +The function :func:`~alchemlyb.visualisation.plot_mbar_omatrix` allows the user +to plot the overlap matrix from +:attr:`~alchemlyb.estimators.MBAR.overlap_maxtrix`. The user can pass +:class:`matplotlib.axes.Axes` into the function to have the overlap maxtrix +drawn on a specific axes. The user could also specify a list of lambda states +to be skipped when labelling the states. + +Please check :ref:`How to plot MBAR overlap matrix ` for +usage. + +API Reference +------------- +.. autofunction:: alchemlyb.visualisation.plot_mbar_omatrix \ No newline at end of file diff --git a/src/alchemlyb/estimators/mbar_.py b/src/alchemlyb/estimators/mbar_.py index 6d3b4e36..171429f6 100644 --- a/src/alchemlyb/estimators/mbar_.py +++ b/src/alchemlyb/estimators/mbar_.py @@ -45,6 +45,9 @@ class MBAR(BaseEstimator): states_ : list Lambda states for which free energy differences were obtained. + overlap_maxtrix: DataFrame + The overlap matrix. + """ def __init__(self, maximum_iterations=10000, relative_tolerance=1.0e-7, @@ -85,6 +88,7 @@ def fit(self, u_nk): verbose=self.verbose) self.states_ = u_nk.columns.values.tolist() + self.overlap_maxtrix = self._mbar.computeOverlap()['matrix'] # set attributes out = self._mbar.getFreeEnergyDifferences(return_theta=True) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py new file mode 100644 index 00000000..ec3771af --- /dev/null +++ b/src/alchemlyb/tests/test_visualisation.py @@ -0,0 +1,17 @@ + +def test_plot_mbar_omatrix(): + '''Just test if the plot runs''' + import pandas as pd + from alchemtest.gmx import load_benzene + from alchemlyb.parsing.gmx import extract_u_nk + from alchemlyb.estimators import MBAR + + bz = load_benzene().data + u_nk_coul = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']]) + mbar_coul = MBAR() + mbar_coul.fit(u_nk_coul) + + from alchemlyb.visualisation.mbar_martix import plot_mbar_omatrix + plot_mbar_omatrix(mbar_coul.overlap_maxtrix) + plot_mbar_omatrix(mbar_coul.overlap_maxtrix, [1,]) + diff --git a/src/alchemlyb/visualisation/__init__.py b/src/alchemlyb/visualisation/__init__.py new file mode 100644 index 00000000..f8cc4c2a --- /dev/null +++ b/src/alchemlyb/visualisation/__init__.py @@ -0,0 +1 @@ +from .mbar_martix import plot_mbar_omatrix \ No newline at end of file diff --git a/src/alchemlyb/visualisation/mbar_martix.py b/src/alchemlyb/visualisation/mbar_martix.py new file mode 100644 index 00000000..6676c739 --- /dev/null +++ b/src/alchemlyb/visualisation/mbar_martix.py @@ -0,0 +1,78 @@ +"""Functions for Plotting the overlay matrix for the Mbar estimator. + +""" +import matplotlib.pyplot as plt +import numpy as np +def plot_mbar_omatrix(matrix, skip_lambda_index=[], ax=None): + '''Plot the MBar overlap matrix. + + Parameters + ---------- + matrix : DataFrame + DataFrame of the overlap matrix obtained from + :attr:`~alchemlyb.estimators.MABR.overlap_maxtrix` + skip_lambda_index : List + list of lambda indices to be omitted from plotting process. + Default: []. + ax : matplotlib.axes.Axes + Matplotlib axes object where the plot will be drawn on. If ax=None, + a new axes will be generated. + + Returns + ------- + matplotlib.axes.Axes + An axes with the overlap matrix drawn. + ''' + # Compute the size of the figure, if ax is not given. + max_prob = matrix.max() + size = len(matrix) + if ax is None: + fig, ax = plt.subplots(figsize=(size / 2, size / 2)) + ax.set_xticks([]) + ax.set_yticks([]) + ax.axis('off') + for i in range(size): + if i != 0: + ax.axvline(x=i, ls='-', lw=0.5, color='k', alpha=0.25) + ax.axhline(y=i, ls='-', lw=0.5, color='k', alpha=0.25) + for j in range(size): + if matrix[j, i] < 0.005: + ii = '' + elif matrix[j, i] > 0.995: + ii = '1.00' + else: + ii = ("{:.2f}".format(matrix[j, i])[1:]) + alf = matrix[j, i] / max_prob + ax.fill_between([i, i + 1], [size - j, size - j], [size - (j + 1), size - (j + 1)], color='k', alpha=alf) + ax.annotate(ii, xy=(i, j), xytext=(i + 0.5, size - (j + 0.5)), size=8, textcoords='data', va='center', + ha='center', color=('k' if alf < 0.5 else 'w')) + + if skip_lambda_index: + ks = [int(l) for l in skip_lambda_index] + ks = np.delete(np.arange(size + len(ks)), ks) + else: + ks = range(size) + for i in range(size): + ax.annotate(ks[i], xy=(i + 0.5, 1), xytext=(i + 0.5, size + 0.5), size=10, textcoords=('data', 'data'), + va='center', ha='center', color='k') + ax.annotate(ks[i], xy=(-0.5, size - (size - 0.5)), xytext=(-0.5, size - (i + 0.5)), size=10, textcoords=('data', 'data'), + va='center', ha='center', color='k') + ax.annotate('$\lambda$', xy=(-0.5, size - (size - 0.5)), xytext=(-0.5, size + 0.5), size=10, textcoords=('data', 'data'), + va='center', ha='center', color='k') + ax.plot([0, size], [0, 0], 'k-', lw=4.0, solid_capstyle='butt') + ax.plot([size, size], [0, size], 'k-', lw=4.0, solid_capstyle='butt') + ax.plot([0, 0], [0, size], 'k-', lw=2.0, solid_capstyle='butt') + ax.plot([0, size], [size, size], 'k-', lw=2.0, solid_capstyle='butt') + + cx = np.repeat(range(size + 1), 2) + cy = sorted(np.repeat(range(size + 1), 2), reverse=True) + ax.plot(cx[2:-1], cy[1:-2], 'k-', lw=2.0) + ax.plot(np.array(cx[2:-3]) + 1, cy[1:-4], 'k-', lw=2.0) + ax.plot(cx[1:-2], np.array(cy[:-3]) - 1, 'k-', lw=2.0) + ax.plot(cx[1:-4], np.array(cy[:-5]) - 2, 'k-', lw=2.0) + + ax.set_xlim(-1, size) + ax.set_ylim(0, size + 1) + return ax + + From b99bda90d8ab13b637ca736a9a4a1eb4b99f8097 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Mon, 24 Aug 2020 16:45:00 +0100 Subject: [PATCH 02/25] depedency --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 0e3e3531..eb689366 100755 --- a/setup.py +++ b/setup.py @@ -29,5 +29,5 @@ license='BSD', long_description=open('README.rst').read(), tests_require = ['pytest', 'alchemtest'], - install_requires=['numpy', 'pandas>=0.23.0', 'pymbar>=3.0.5', 'scipy', 'scikit-learn'] + install_requires=['numpy', 'pandas>=0.23.0', 'pymbar>=3.0.5', 'scipy', 'scikit-learn', 'matplotlib'] ) From 6111a8ede979cef02b9e32aaf1605f4204a98efc Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Mon, 24 Aug 2020 16:53:12 +0100 Subject: [PATCH 03/25] py2 compatibility --- src/alchemlyb/visualisation/mbar_martix.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/alchemlyb/visualisation/mbar_martix.py b/src/alchemlyb/visualisation/mbar_martix.py index 6676c739..889902ed 100644 --- a/src/alchemlyb/visualisation/mbar_martix.py +++ b/src/alchemlyb/visualisation/mbar_martix.py @@ -1,6 +1,7 @@ """Functions for Plotting the overlay matrix for the Mbar estimator. """ +from __future__ import division import matplotlib.pyplot as plt import numpy as np def plot_mbar_omatrix(matrix, skip_lambda_index=[], ax=None): From bf80720420af898fbea4aa095633d104699c60c4 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Mon, 24 Aug 2020 16:59:52 +0100 Subject: [PATCH 04/25] backend error --- src/alchemlyb/visualisation/mbar_martix.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/alchemlyb/visualisation/mbar_martix.py b/src/alchemlyb/visualisation/mbar_martix.py index 889902ed..e854d74f 100644 --- a/src/alchemlyb/visualisation/mbar_martix.py +++ b/src/alchemlyb/visualisation/mbar_martix.py @@ -2,6 +2,8 @@ """ from __future__ import division +import matplotlib +matplotlib.use('Agg') import matplotlib.pyplot as plt import numpy as np def plot_mbar_omatrix(matrix, skip_lambda_index=[], ax=None): From c5ac34a3cf1ab1474467cf2b3bae7e27cd8e0f5c Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Mon, 24 Aug 2020 17:06:26 +0100 Subject: [PATCH 05/25] bump coverage --- src/alchemlyb/tests/test_visualisation.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index ec3771af..59031681 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -15,3 +15,9 @@ def test_plot_mbar_omatrix(): plot_mbar_omatrix(mbar_coul.overlap_maxtrix) plot_mbar_omatrix(mbar_coul.overlap_maxtrix, [1,]) + # Bump up coverage + overlap_maxtrix = mbar_coul.overlap_maxtrix + overlap_maxtrix[0,0] = 0.0025 + overlap_maxtrix[-1, -1] = 0.9975 + plot_mbar_omatrix(overlap_maxtrix) + From 8e554095d4283de45db15f82a1945f73da954fb3 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Mon, 24 Aug 2020 18:13:58 +0100 Subject: [PATCH 06/25] citation --- docs/visualisation.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/visualisation.rst b/docs/visualisation.rst index c84cfcae..f2afc938 100644 --- a/docs/visualisation.rst +++ b/docs/visualisation.rst @@ -11,7 +11,7 @@ The accuracy of the :class:`~alchemlyb.estimators.MBAR` estimator depends on the overlap between different lambda states. The overlap matrix from the :class:`~alchemlyb.estimators.MBAR` estimator could be plotted to check the degree of overlap. It is recommended that there should be at least -**0.05** overlap between neighboring states. :: +**0.03** [Klimovich2015]_ overlap between neighboring states. :: >>> import pandas as pd >>> from alchemtest.gmx import load_benzene @@ -26,3 +26,7 @@ the degree of overlap. It is recommended that there should be at least >>> from alchemlyb.visualisation.mbar_martix import plot_mbar_omatrix >>> ax = plot_mbar_omatrix(mbar_coul.overlap_maxtrix) >>> ax.figure.savefig('O_MBAR.pdf', bbox_inches='tight', pad_inches=0.0) + +.. [Klimovich2015] Klimovich, P.V., Shirts, M.R. & Mobley, D.L. Guidelines for + the analysis of free energy calculations. J Comput Aided Mol Des 29, 397–411 + (2015). https://doi.org/10.1007/s10822-015-9840-9 From 5f03a3e44f4a0bfb9ced121543f300b5098fa1c8 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Tue, 25 Aug 2020 20:56:35 +0100 Subject: [PATCH 07/25] apply change --- .gitignore | 2 ++ .travis.yml | 4 ++++ docs/api_principles.rst | 2 +- docs/images/O_MBAR.png | Bin 0 -> 9960 bytes docs/index.rst | 1 + docs/visualisation.rst | 11 ++++++--- ...isualisation.plot_mbar_overlap_matrix.rst} | 12 ++++++---- src/alchemlyb/estimators/mbar_.py | 9 +++++-- src/alchemlyb/tests/test_visualisation.py | 22 +++++++++++------- src/alchemlyb/visualisation/__init__.py | 2 +- .../{mbar_martix.py => mbar_matrix.py} | 8 +++---- 11 files changed, 48 insertions(+), 25 deletions(-) create mode 100644 docs/images/O_MBAR.png rename docs/visualisation/{alchemlyb.visualisation.plot_mbar_omatrix.rst => alchemlyb.visualisation.plot_mbar_overlap_matrix.rst} (52%) rename src/alchemlyb/visualisation/{mbar_martix.py => mbar_matrix.py} (95%) diff --git a/.gitignore b/.gitignore index 2a45ed72..f71e81e2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ *.pyc .vscode +*.DS_Store +build diff --git a/.travis.yml b/.travis.yml index 3d9c07ef..d4f03437 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,6 +2,10 @@ language: python sudo: required dist: xenial +env: + global: + - MPLBACKEND=agg + python: - "2.7" - "3.5" diff --git a/docs/api_principles.rst b/docs/api_principles.rst index e644fb37..e0c6d73c 100644 --- a/docs/api_principles.rst +++ b/docs/api_principles.rst @@ -1,5 +1,5 @@ API principles -============ +============== The following is an overview over the guiding principles and ideas that underpin the API of alchemlyb. diff --git a/docs/images/O_MBAR.png b/docs/images/O_MBAR.png new file mode 100644 index 0000000000000000000000000000000000000000..65effe53df2696df1c618d0fcda6e4877756c61f GIT binary patch literal 9960 zcma)ibyU>-w)fB?rJ}S9JcQCINP{SiC=DVZU4nF{iU^2+lt?2nG>Fn6B`w_{A>Gn- z_dMsWv+jA{bMCr-80gGzzOnbG_5?qAC`)vi>M{a>Ad;86uL^%Qz_0pCxbS-oEuIkk zA>{Nx+eyvN%*oZr!4#ov!O_z0IUgq))U?JG@Q>blJPjn@7EI< zU*&1NA$a?${n2F{%HIK&Ll#u&vnf)`^Zd;t+QOfI^QLK9h0gN6jShZN!FD(EWryMW zfVpdbvA(V-(&wDmJqw%`Gibx9pxG#!J`+Wf(2{HjE&AvqIS-tNdKS+jO>|<49lDdm z6kFfjd{rC0GOn$3#B(pJAK#m8h@AO~w%B)U%v{+xNyEKMq?ht&jWGQO8I;cZu4#$kG+*wrKNl-qTASwuZb9L-DrTTE&XX9=b}Ryn48 zFT5Wni3S%Ia`i}Dhz=JT2!{B`%j4J9)_NT;XPNb$+{wwwX=-YE9~09uRpZ&ax<^RC z?H3t&t#rxL-JQpEL+|0IyA9iOt?v7y=D}~?%(NnX`9sbpX+&H?Z?I|pOqV6;r+U9C zWoemXV{FV|H&tCB-W7LWQqqrxlr+7%0ELPa6BAnya{N6|k&pB_Ox5?sn4h0FiJCSr z@Yac0SRF3x@cZ=fqt#XPPi9>o+b&KQCAv@s>gtTHOFyNis$2_cXlU3~;(0%J+YYRL zOerb3CMzpjQc|+guauaRojubOO#1TGE0zu$8=Ff61anKh>1OL4bKlp;zhPLwaz@K? zQK*sf7Ys~HOn){e$MeQ`cz6h@__JGFQwuCs`uz{4Ji}*apOq%06cxQkp;#v-Cal<9 z9UU1&M8@;H-D$Wf2$d3rqielS`tqL$mM)u_ne|n<+O>Cdyl`;%OGqf5^J!(AvhQ3= z*e6lKYAJD+1pDL;Gm>_>b1?AEVFv%&R&98yY3%86R< z1my&QtCW=TR#qsT(D$*itshv_Myg!7h>3~!tU5CK=>r=iot*fKi;JuG$E;5e>MkA> z4iD!tI5!Lq-gxrl$>qzJt=PG4->#pWoU}P}qj?h+W(w6yBkFdAi;Jt&%DPNaNrho5GS_f|#E{`6b$Q}I&~K@qU)lwgO+NNZ>eGJG?Q#E;>8 zs$gV9UuxEs=6SHOwC7Qp^xC3K*JFRJ?fmQ{+7vg?|Lo)-;q23Z>W>m78yohvwl3yX`p-|Ei7 zGBOyc1fMGf5#KEO_U#rMTgb*_)o__ml>2J_w>P1oP4HGfKg8_Z919DJ(M0(RvtBMv z&Vm=#hqKKgCMG67tKA)|Yig+VB$%0B6HyBkIj?GdN=&S9S*JoMf8;-rvv6&VB>0dLJ429X^LzkGg(z^eAaHxh0IcKkY$Cr2+Zio=YFx(?YC)Pelu_ zgoMOVq3@Y>qMh}SQ;_Nl3<63yIbEKJP#FJJhh#&{OUCL03rOlOjkfBcXuU7A** zeosnD+Lb2ZbLv(B7mR-Oldm=`L9bu8d+x90IErMhOp~bDMGP&=6&be@6Vu$ShwHIo zhX${D!1>PD*f=;W>`oP@q@3KP5}YjaC#uwj3+>TlgLOza0_K{+hs+eBZrp;#wQ2vh96*v);MlO5iMF&^iBMpsL`S~0L1qEiqTT?Hg5!>hekWn{eqE90Z1I!bl=j-jg^RRDDL{)V%%r*W&YksNogb#!%g z*$mRt;y`njTF1o1VA9E7>-u_Z zs07iHLyJ>L_*9QNYZd6p>*!Fzy&0Vy@0ktPIj@qzo|WRb{mD3VVx-K~@*{8BaZAjwcsOP?n%rj$S>4y)Ap}EyicXFkmwB+TMxpK%4)5*_8$=}}V>Z%2+72OgvoNo$Xenv7Oz;&wVBFYo~L_3&Yu_ z;81G8cK}R9UPlfzH*a?FQf#RB9+NKiBo`L5v$Kc3m}}_;5UX@qZx1HDb-0usu0F`Y z!NJ7C6ERw9&Sy6zI`-98QBo3t$k#2K^*vwo)h;r`K>D8YK$||B{}F*&5fK*N8M7|j z?iSe;wNC&{7Znv%AFN5Jt`>!>Y9crAYVgOxD}OPEOD*^uJ0&kKbZN-~O-lP~qmA8( zBHITCS{f39g0b*{yX)Vk9d_pB3k>SIR|a!2mWJ}F@rY?=rZ3I~r@aquzIyd4TpiiQ zNUF_?Oi_&GB%~Gy1*#j;-dJ7j_|+K1!pVs{Gj3OBx;9dbvNEXkiigic40?V%Jw2sz zM|}PIwOzlpw>RR+x92%1J9e?d6xPm)iV7Mr4|V6=$u3w8J`vGp(~cOR+3~S4l=yYy)uH^CK|yoBf4^pwN1twC-8ohEd1shx#@n}-fenzBqLK2< zZrb90#mr{XDZLhQa&mOS!pBLdMj17#lRqoDg-`)~hJk&6PTRY?YOI@pfR0-;QsRd* zSU`Obpxq6;_cQ`6ro-FY709^sJE`reA7?%L;tbu-_Qh&Y!F_#Up*fmO3xO^q-TavB z?Ay2`w4Glp`#&VU*42Z359L2$E3B!hsqi`~**Wg}toUbjIDCJ7d~SC3`%8RME8AT` zXV;~(gK6{K-(3g!eu;pC0XPIB#m0D+mX@bba*}(o(zei$Z{EE#o38bqUs%BSnf8Fd z%iB9%$dQ@Ptdow3soAF4A7iM{KqXn+yHl0$P)c1>b7#_Zdiz%pvc&tuWnywNp2zsL z&dJRBdKcVJTX(liRzHhuPLNR~K6^CD6&f1!-Wr-i$c4NPxj8vG^{U+-#}k&?jPimI zdhqDc73g~S=*Q~n8CJ-(_4RmBci#B;cwn|iI_OJ>LnwD#Bu2oZzgQ0!NFB7hM#aas z!A^riLruoMSbm6!xkgEubIcN-m38ar=%{_QBbIA!c{wMK;PU0~lan{)Bj3DvL&v}% zt>Fz^Y1)c4=y`s!z*E7RqGiV!yM{#22LeH2D5csRM~lQki<0;PGkU`larB;A3p|FprEKo$!A8Gnwn}RwEi|E zgn^&`0{~*CQv7$I6+%fvL%MXCFicc9cGN8uGr*0Jk&&6-zca^fRtipgMQv@_>%0BJ znesk(`n|CcV`;cBdFKH3_y)8knt?PlZh*>I1av4UTleTys+bK55H}1vt&Q~NeRnmsw$AN5_4rvm`Z}|^TJLCb;sX0n zVopBtWmHtu=w3)jNQ5b+eB>o!Vu7mjL=o59Y;1EmSxLN};(;#$Y9FWs={MpU8g9kd~g_RZ#b=`7&(9VxpV}Ez3vurad>s-oAbNo>6`e|895ydj}{( zKmm9`u})c7D&y>T?f+_h zW5GwT4c2w~f8NCaS=uwx5E$LwpQVy)2D}O?GV5WYRcH;U>NL=3&>*Nx-~Bb6uXfY# z9%e~LGFCb*e}`i(DlV4F>Myn$WZl`>0danyzOES(5P*p$9e9(yqoeV1-{%7xh7MkM z_98H_*LE{G>VO`kdQ;PV6P--0Lj9)ZW+QX+paw%aaq%P&6lehiCJEl4#i}83DN2dk zent|Pn8uif}dHw$LsN?WOALB|#Sh z2{rfjJ}ztbBw1+$NM_;SFkK%jJD43PW_F5dCOkZ?^90mdHWsy*qV~we}h&8{sD4{<6CM1O5rm)k@;r1MA z1za9WbI|(Q+SbmFs?4kvKb5e{+Q`gYGLY2rSlKF8zk!3Jqr9S`by4U`@00zjG&G^P zxonY)+kJ{$#n!`|DaUe@6cld~6K_T`_I|T7s6G0fpgVr596$4%J&|N)Z7mee74Ep! z{3#qAz~f^~X{rDb*NyiP5!-S&=3Ng=O__kwnFIthDPz+G;BxXKn$SRP^6c>N@cqY+ z3NK#pc+@0`dG6mH#V2V@!&*}M{ME}39`Uoz*B5lVte02y$}!mP`wt$xn5p-tGW{Ng zjg9^O)2F=3f{2|R2l0zj$9vM!93N%3tRo{Mm-@3v@}GSB`lDzK%(%neibh&m8rpks zbHD%c<%{n0K?0xD(;W<%cQiTbw>8xygoQum<*}=(scq^^7d#bNTU*4{&lG?!lSnloZH8D3Yb={mgc@E%v zAH?6@?fo-JC*hkwX@JPGd-vt;l#7r6EBFDbs_oIq8nYb3qtOK-U^|xKMJAj`ocJ&Z zE*`Ll5k%gwfp5)_fHJStix)4#M3CfM`iu+=flx?-KYoq&9U%&QGp&R$*j9g5rY*I1;B@SO^~<9|#f< zejw6RRaNh~y9ex6rKE8JIX1)m{NbDUryR9IJQ>hz?cs!23l>!0- z@$vDEpo0{crGSYjd7fcEI;;6~_r+h4jPxD6z$YFV!S|RKNHa1rCJr9a0Lym{xnBqQ z*_83|LNKSLtxbOK$Ea+>-JDF2(Wk2h7ra&Hu(i<|&uYZ(*kJ?SoowU|5wBJB)TXL@hKWk& zj-4ZH+`xusp7hxd-|%eUcRi(|$>zc7zjh1}27>J5Buw@%YHI4z1Ufdh+%ar3jk6mZ zPrlBBMC7*~3f)^Bev&}~?DQLM6U2-TBhv8r&+>J0^8WG{##8Q#H%nK2Pq$k?eE9Ga z{=uv>&dB8qRLo;Nz0k0*F7FUBenu)~@G*K61phAN`lq#OWhCv4=L^!W@%Xj2)=8=& z+E8TFXr*pKcF9LE?oSc>ojcKSaTI7Z2Mz9fwrW7hZHNFMc>Xi%9S9B2 z8{J^JJbDC5*h#^XLcxrDv5Z0;$8_VB+yBCZLW54g%r z#fvKS92|IQX=zdBSQL94@gzTiRRJvIK~e*zqvai3_4a_83d{H>ac`<8PsVY31^12q zEcJe|=*6tjSGji&^Oi~yBpb*unMOFat*xzkLvHvBE-|fWw)zh>EJzIXy|y0!MQ^Zp zfF3mhGuISK$yZcXhNgtI@o$Ag{_ZF}Qb>`gQyP$-P>TJKPXT3bREdMDLhBQ4TQ5L; zXyV(+#Q@M78ygeV?rJ&)dU~Azd{j>SX1SXA7QJb&`tX3W7x*kjOFpX#Ou8|aSNU4a z%+KEkv-jakT2)&|M`LR%UYJZEMFKER0|^#!8=PsCFtB|?BB8gZ}kTIInIsBw5C5}o=~3)F&|0Q;p= z=Ocz*h(@j^0|t#!9i6ti^F1ygIg6p2#^7%%z@Wmt zN#_oN%We#pIM-f1SuqV58PNtFY=C0|VnnH^?X_rA`=0EtE7+a_d4o4o0t2qmU+g6A z2ikbD(ivTVkP)(WaOqe7fZVqUlwY)IviCkEoK``^+TV1Mr&keK@U&{mkw{xZqZuMN zfM27fpQ)&oA#j!V@87?2^(uy>heu^elp6(eJZQ@)w>hH7s3<_$R1or+{OPQ&cz37rKHe-sRm3_wVDg33V>2xzJBenGQh&r+CUMt2iksTHiT>FiSORZ;Nj8H zFSrlS&6*%sB7_IurhO9Nvti{BHnRY7w1S<1cvM@P0?4N9t8E|69G;LU4^Qz!BO-KY zfyP6-sfpPrE08kz9;#uWnIGI5gfM+8-1oFGZ>-e3=lAzl1ewXi_}CY9bS=dA|0~}8 zYw~a0>xDTtk_xtcd(=ExCxZk|Xcn#jlC20YKgg3!Dc&;Wss@_wUQO}%ogxHM?lTH? z8SMv(0sHj}@lfXUpEEMFGe}x|{n$IM-?#w}k=Z!EryEe!(ewmU(fgdN>lSx6sk4)l z30QZsWoZzBIwdCfk&N-)ry`>IN0(wAuDJi4*nSN|3i^%`(_Q>fuvfREzIAnhPQHEn zHr?I3+Jlr9;q7re#^uhdG0nReRClbJkQaxtKmzASiyxw+f6dIKv`>UoqGMBtUpP5K zaWMb?(Z>HoeYEDG74v|RL;whi8qWhWu%=cA%5>ZWAoX{5cAkMrc3S$0p%VkVLUi@& zRkW$J+``N=@Ht`v$s!wL8W0c=cgH$L^+NpCtsn?LRkgJ*)kBP-o_}5I(W80@V$tCL z{P}bCrOBb{YrHBhVz|?(6I}IB!PtVSB zP^ejen7?x|UY9j(aK<{|)V;p}7ZCvR!s7Kc{$6;>>H0tH#Mbd@VSla`S6_cW1SHHN zB2PwLAi6`x{2EUs=;{4|vml(etJx!yl9Gb3>N*mMgmY;8B;*K9U%mC603C#Z3<~`S9T%61 zP@k5jCTsAWN8Ya;HYS8jB`ttYVWkj1l;HF7@qv^3J~fpDLNwyjC#6`rDQWJp_BjgtFSp6aXqbk86HjoF_Hoj3jAey`keGQ&r#hVkPc{t z!}F6(bbb^O5dm~OHt>5z86*ihR|M-zKmN_G6~-rdSbhJn20dbB$J7C>^wjmRb`Ss# zGxZqYG&gjRkQt*7XITW225RaK(EV*)T{x&h{DZYoX&8Y3n7_p!!fliu=XJlH1ELxD z2=QlaGzx}7#ctcrH1o8LtfI*vzL#Z@2pwnZ{w z8&kB7g_bU08lFd1}4_ zqX+mQz2%`4QOmoV9Uws@(F$n%D+jy`Jr0v}R*&wz0_yESBqJDd@V#@8Qv8Em@JX4A zd@pt(YSo3!(MH>=OUBa%Wt7~=DHu0e>RIlsCgys9k;blex z>@XD@%F{6$7K{nxf&d8or(jj>i*q|LU65Jx<$irJuxlKpVwi*B)ZE;B3M~V$pV}Me zn+&N3j(S1^tglR2&d$P$`e6EzHh*%xHhx1xgH$}Z;d2EA1xU1h{T(WKdm*&6w9sA} zuA`HeV|N(x4B4Z}YrYqg23|0YEVZI`l?-oQh775!8_V`KAud)J)txjT0W$oyOC%(V zYon!LbS7NCdr0qQNt^{k5rhz$;%a|;eiOX($k-U`E@|fjn1X!%{2B5ly8q3Mt=I;D z`2S1(8eEn}stBh4QOcS8%Gp`jR9IFIMP z4M=9>SZ14S0Eh;Ci2tu$yJj(z7oXEL5IixVXV#O{1s&A^jZoopCJZ(kHwCmdderKB zaWd_7FzI@oj0~lZWRZ|a?h!u;1^ruWI}SWh-T(&Be?>1~ZOsOdZKT?rPe1^Au1HG5 z`@h`<^TSVWJ(Z0?B#z5{OyC3?VBX@&C|{ZNFpG*}7Xv+)5I^0LhHxV5gM{F3tjP{JwKH;RjR5YH%Jbzi>tVl@ z8biQf=rxeRFziY=``t6?-F`KnK#9MJ=<)=gcD{RM?9GCeh%C;^2OrgcBp5856%I zzNr8c8Fb*l%K8>QGR|4u1UwD|d2E2?#^C{fdC{>N3QLG<)WP*97K<(PZLYtRY4&H(DLKpp$_3VjYJuh}RSafo5>j7Us( zbcBSEbFmhpSvWR!SZo|7sxU&auRUI21G%31u{J13Hcav@n^wVX_e+tCJeGZoKo8jU z;DjuO3lgoy%Af*(Ccyt7-X|sLGYrZq4|L;kN0*fFx@=B$5pY%4L;Q!(E`IhZGJ?mj z9>J+s{>WOyc}fG+O=Nt$yo}82ow?1;p4`G(1sNHzo@a0r?eKgM!C<1Uk&;r1d)-Bl zbLwUD$S=Z_C zAv*$OHZn2-kCl@~LX<$Vb^CPP&X?iNof`n~J^(Am+!PQrFacb?e2JZMOh`vhj}h3A zGLb0awDc(7HsOQHmi4`R_qJe)`JJBjhOxMRP*70wPB*M-6q*O*=hq(@A0d>FVb}SK zc~IM4HH=k7btCo()0I=rE50LCcs*+K$15&$|l zVw(uv6xbRF)?d&t5drlr5i!u`pm9rAFJOMz?uBrSon`0ZYBl5*rDeWyi9Pu4&1>3X b7X>> mbar_coul = MBAR() >>> mbar_coul.fit(u_nk_coul) - >>> from alchemlyb.visualisation.mbar_martix import plot_mbar_omatrix - >>> ax = plot_mbar_omatrix(mbar_coul.overlap_maxtrix) + >>> from alchemlyb.visualisation import plot_mbar_overlap_matrix + >>> ax = plot_mbar_overlap_matrix(mbar_coul.overlap_matrix) >>> ax.figure.savefig('O_MBAR.pdf', bbox_inches='tight', pad_inches=0.0) + +Will give a plot looks like this .. image:: images/O_MBAR.png + .. [Klimovich2015] Klimovich, P.V., Shirts, M.R. & Mobley, D.L. Guidelines for the analysis of free energy calculations. J Comput Aided Mol Des 29, 397–411 (2015). https://doi.org/10.1007/s10822-015-9840-9 diff --git a/docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst b/docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst similarity index 52% rename from docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst rename to docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst index ccdaa6af..5323d1a9 100644 --- a/docs/visualisation/alchemlyb.visualisation.plot_mbar_omatrix.rst +++ b/docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst @@ -1,9 +1,11 @@ -Plot Overlap Matrix from MBar +.. _plot_overlap: + +Plot Overlap Matrix from MBAR ============================= -The function :func:`~alchemlyb.visualisation.plot_mbar_omatrix` allows the user -to plot the overlap matrix from -:attr:`~alchemlyb.estimators.MBAR.overlap_maxtrix`. The user can pass +The function :func:`~alchemlyb.visualisation.plot_mbar_overlap_matrix` allows +the user to plot the overlap matrix from +:attr:`~alchemlyb.estimators.MBAR.overlap_matrix`. The user can pass :class:`matplotlib.axes.Axes` into the function to have the overlap maxtrix drawn on a specific axes. The user could also specify a list of lambda states to be skipped when labelling the states. @@ -13,4 +15,4 @@ usage. API Reference ------------- -.. autofunction:: alchemlyb.visualisation.plot_mbar_omatrix \ No newline at end of file +.. autofunction:: alchemlyb.visualisation.plot_mbar_overlap_matrix \ No newline at end of file diff --git a/src/alchemlyb/estimators/mbar_.py b/src/alchemlyb/estimators/mbar_.py index 171429f6..9a2fb519 100644 --- a/src/alchemlyb/estimators/mbar_.py +++ b/src/alchemlyb/estimators/mbar_.py @@ -45,7 +45,7 @@ class MBAR(BaseEstimator): states_ : list Lambda states for which free energy differences were obtained. - overlap_maxtrix: DataFrame + overlap_matrix: DataFrame The overlap matrix. """ @@ -88,7 +88,6 @@ def fit(self, u_nk): verbose=self.verbose) self.states_ = u_nk.columns.values.tolist() - self.overlap_maxtrix = self._mbar.computeOverlap()['matrix'] # set attributes out = self._mbar.getFreeEnergyDifferences(return_theta=True) @@ -102,3 +101,9 @@ def fit(self, u_nk): def predict(self, u_ln): pass + + def _get_overlap_matrix(self): + return self._mbar.computeOverlap()['matrix'] + + overlap_matrix = property(_get_overlap_matrix) + diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index 59031681..5de33ab5 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -1,23 +1,27 @@ +import matplotlib +import pandas as pd + +from alchemtest.gmx import load_benzene +from alchemlyb.parsing.gmx import extract_u_nk +from alchemlyb.estimators import MBAR +from alchemlyb.visualisation.mbar_matrix import plot_mbar_overlap_matrix def test_plot_mbar_omatrix(): '''Just test if the plot runs''' - import pandas as pd - from alchemtest.gmx import load_benzene - from alchemlyb.parsing.gmx import extract_u_nk - from alchemlyb.estimators import MBAR - bz = load_benzene().data u_nk_coul = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']]) mbar_coul = MBAR() mbar_coul.fit(u_nk_coul) - from alchemlyb.visualisation.mbar_martix import plot_mbar_omatrix - plot_mbar_omatrix(mbar_coul.overlap_maxtrix) - plot_mbar_omatrix(mbar_coul.overlap_maxtrix, [1,]) + isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_maxtrix), + matplotlib.axes.Axes) + isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_maxtrix, [1,]), + matplotlib.axes.Axes) # Bump up coverage overlap_maxtrix = mbar_coul.overlap_maxtrix overlap_maxtrix[0,0] = 0.0025 overlap_maxtrix[-1, -1] = 0.9975 - plot_mbar_omatrix(overlap_maxtrix) + isinstance(overlap_maxtrix, + matplotlib.axes.Axes) diff --git a/src/alchemlyb/visualisation/__init__.py b/src/alchemlyb/visualisation/__init__.py index f8cc4c2a..4fc73051 100644 --- a/src/alchemlyb/visualisation/__init__.py +++ b/src/alchemlyb/visualisation/__init__.py @@ -1 +1 @@ -from .mbar_martix import plot_mbar_omatrix \ No newline at end of file +from .mbar_matrix import plot_mbar_overlap_matrix \ No newline at end of file diff --git a/src/alchemlyb/visualisation/mbar_martix.py b/src/alchemlyb/visualisation/mbar_matrix.py similarity index 95% rename from src/alchemlyb/visualisation/mbar_martix.py rename to src/alchemlyb/visualisation/mbar_matrix.py index e854d74f..14fb3210 100644 --- a/src/alchemlyb/visualisation/mbar_martix.py +++ b/src/alchemlyb/visualisation/mbar_matrix.py @@ -2,18 +2,18 @@ """ from __future__ import division -import matplotlib -matplotlib.use('Agg') + import matplotlib.pyplot as plt import numpy as np -def plot_mbar_omatrix(matrix, skip_lambda_index=[], ax=None): + +def plot_mbar_overlap_matrix(matrix, skip_lambda_index=[], ax=None): '''Plot the MBar overlap matrix. Parameters ---------- matrix : DataFrame DataFrame of the overlap matrix obtained from - :attr:`~alchemlyb.estimators.MABR.overlap_maxtrix` + :attr:`~alchemlyb.estimators.MBAR.overlap_matrix` skip_lambda_index : List list of lambda indices to be omitted from plotting process. Default: []. From c327215f9371846be7df23aa5f367350d5971079 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Tue, 25 Aug 2020 21:03:17 +0100 Subject: [PATCH 08/25] fix test --- src/alchemlyb/tests/test_visualisation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index 5de33ab5..562609ed 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -13,13 +13,13 @@ def test_plot_mbar_omatrix(): mbar_coul = MBAR() mbar_coul.fit(u_nk_coul) - isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_maxtrix), + isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_matrix), matplotlib.axes.Axes) - isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_maxtrix, [1,]), + isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_matrix, [1,]), matplotlib.axes.Axes) # Bump up coverage - overlap_maxtrix = mbar_coul.overlap_maxtrix + overlap_maxtrix = mbar_coul.overlap_matrix overlap_maxtrix[0,0] = 0.0025 overlap_maxtrix[-1, -1] = 0.9975 isinstance(overlap_maxtrix, From 095868e111404fe11241d5c358660dc7e135f164 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Tue, 25 Aug 2020 21:45:59 +0100 Subject: [PATCH 09/25] fix test --- src/alchemlyb/tests/test_visualisation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index 562609ed..ff60b6f4 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -22,6 +22,6 @@ def test_plot_mbar_omatrix(): overlap_maxtrix = mbar_coul.overlap_matrix overlap_maxtrix[0,0] = 0.0025 overlap_maxtrix[-1, -1] = 0.9975 - isinstance(overlap_maxtrix, + isinstance(plot_mbar_overlap_matrix(overlap_maxtrix), matplotlib.axes.Axes) From 0b998b910265335ff4ff592924d5bd02a4f5e598 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Tue, 25 Aug 2020 22:01:09 +0100 Subject: [PATCH 10/25] fix test --- src/alchemlyb/tests/test_visualisation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index ff60b6f4..d482caad 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -13,15 +13,15 @@ def test_plot_mbar_omatrix(): mbar_coul = MBAR() mbar_coul.fit(u_nk_coul) - isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_matrix), + assert isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_matrix), matplotlib.axes.Axes) - isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_matrix, [1,]), + assert isinstance(plot_mbar_overlap_matrix(mbar_coul.overlap_matrix, [1,]), matplotlib.axes.Axes) # Bump up coverage overlap_maxtrix = mbar_coul.overlap_matrix overlap_maxtrix[0,0] = 0.0025 overlap_maxtrix[-1, -1] = 0.9975 - isinstance(plot_mbar_overlap_matrix(overlap_maxtrix), + assert isinstance(plot_mbar_overlap_matrix(overlap_maxtrix), matplotlib.axes.Axes) From a98027f8625e53b33621c9f37ee70688797c1607 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Wed, 26 Aug 2020 10:31:21 +0100 Subject: [PATCH 11/25] fix visulisatiuon --- AUTHORS | 4 +++- docs/conf.py | 6 ++++-- docs/visualisation.rst | 4 +++- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/AUTHORS b/AUTHORS index 672262b7..c48f9fa2 100644 --- a/AUTHORS +++ b/AUTHORS @@ -33,4 +33,6 @@ Chronological list of authors 2019 - Victoria Lim (@vlim) - Hyungro Lee (@lee212) - - Mohammad S. Barhaghi (@msoroush) \ No newline at end of file + - Mohammad S. Barhaghi (@msoroush) +2020 + - Zhiyi Wu (@xiki-tempula) \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py index 517695f4..28d24289 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -50,8 +50,10 @@ # General information about the project. project = u'alchemlyb' -author = u'David Dotson and contributors' -copyright = u'2017-2019, ' + author +author = u'''David Dotson, Ian Kenney, Oliver Beckstein, Shuai Liu, +Travis Jensen, Bryce Allen, Dominik Wille, Victoria Lim, Hyungro Lee, +Mohammad S. Barhaghi, Zhiyi Wu''' +copyright = u'2017-2020, ' + author # The version info for the project you're documenting, acts as replacement for diff --git a/docs/visualisation.rst b/docs/visualisation.rst index 2a848f14..e7a03336 100644 --- a/docs/visualisation.rst +++ b/docs/visualisation.rst @@ -30,7 +30,9 @@ the degree of overlap. It is recommended that there should be at least >>> ax.figure.savefig('O_MBAR.pdf', bbox_inches='tight', pad_inches=0.0) -Will give a plot looks like this .. image:: images/O_MBAR.png +Will give a plot looks like this + +.. image:: images/O_MBAR.png .. [Klimovich2015] Klimovich, P.V., Shirts, M.R. & Mobley, D.L. Guidelines for the analysis of free energy calculations. J Comput Aided Mol Des 29, 397–411 From 5613966b92b8cff8e02ba12c289e503ffae4797d Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Wed, 26 Aug 2020 10:33:16 +0100 Subject: [PATCH 12/25] Update CHANGES --- CHANGES | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGES b/CHANGES index 98638667..fcc8e133 100644 --- a/CHANGES +++ b/CHANGES @@ -15,11 +15,12 @@ The rules for this file: ------------------------------------------------------------------------------ -??/??/2020 wehs7661, dotsdl +??/??/2020 wehs7661, dotsdl, xiki-tempula * 0.?.? Enhancements + - Allow the overlap matrix of the MBAR estimated to be plotted. (PR #107) Deprecations From 505e0d66843c6c339da3710a6ecbc904c95a0450 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Wed, 26 Aug 2020 10:34:03 +0100 Subject: [PATCH 13/25] Update CHANGES --- CHANGES | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES b/CHANGES index fcc8e133..b95aa851 100644 --- a/CHANGES +++ b/CHANGES @@ -20,7 +20,7 @@ The rules for this file: * 0.?.? Enhancements - - Allow the overlap matrix of the MBAR estimated to be plotted. (PR #107) + - Allow the overlap matrix of the MBAR estimator to be plotted. (PR #107) Deprecations From 453242cb6d688fed9e1b8a1eed25cfb29341a638 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Wed, 26 Aug 2020 19:38:04 +0100 Subject: [PATCH 14/25] Update src/alchemlyb/estimators/mbar_.py Co-authored-by: Oliver Beckstein --- src/alchemlyb/estimators/mbar_.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/alchemlyb/estimators/mbar_.py b/src/alchemlyb/estimators/mbar_.py index 9a2fb519..6aa1821b 100644 --- a/src/alchemlyb/estimators/mbar_.py +++ b/src/alchemlyb/estimators/mbar_.py @@ -102,8 +102,18 @@ def fit(self, u_nk): def predict(self, u_ln): pass - def _get_overlap_matrix(self): + @property + def overlap_matrix(self): + r"""MBAR overlap matrix. + + The estimated state overlap matrix :math:`O_{ij}` is an estimate of the probability + of observing a sample from state :math:`i` in state :math:`j`. + + The :attr:`overlap_matrix` is computed on-the-fly. Assign it to a variable if + you plan to re-use it. + + See Also + --------- + pymbar.mbar.MBAR.computeOverlap + """ return self._mbar.computeOverlap()['matrix'] - - overlap_matrix = property(_get_overlap_matrix) - From 4492fabe8336f4f26f0e9ceb56c3fa4c15b44fdf Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Wed, 26 Aug 2020 19:47:10 +0100 Subject: [PATCH 15/25] update docs --- src/alchemlyb/visualisation/mbar_matrix.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/alchemlyb/visualisation/mbar_matrix.py b/src/alchemlyb/visualisation/mbar_matrix.py index 14fb3210..b5ee8c5d 100644 --- a/src/alchemlyb/visualisation/mbar_matrix.py +++ b/src/alchemlyb/visualisation/mbar_matrix.py @@ -1,4 +1,12 @@ -"""Functions for Plotting the overlay matrix for the Mbar estimator. +"""Functions for Plotting the overlay matrix for the MBAR estimator. + +To assess the quality of the MBAR estimator, the overlap matrix between +the lambda states can be computed and the more overlap is observed between +the states, the more reliable the estimate is. One way of accessing the +quality of the overlap matrix is by plotting it. + +The code for producing the overlap matrix plot is taken from +: `Alchemical Analysis `_. """ from __future__ import division @@ -7,7 +15,7 @@ import numpy as np def plot_mbar_overlap_matrix(matrix, skip_lambda_index=[], ax=None): - '''Plot the MBar overlap matrix. + '''Plot the MBAR overlap matrix. Parameters ---------- @@ -25,6 +33,12 @@ def plot_mbar_overlap_matrix(matrix, skip_lambda_index=[], ax=None): ------- matplotlib.axes.Axes An axes with the overlap matrix drawn. + + Notes + ----- + The code is taken and modified from + : `Alchemical Analysis `_ + ''' # Compute the size of the figure, if ax is not given. max_prob = matrix.max() From 8cfd5204468de8cbe7ae698accaff1521dcf683e Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Wed, 26 Aug 2020 20:26:54 +0100 Subject: [PATCH 16/25] Change data type. --- src/alchemlyb/estimators/mbar_.py | 2 +- src/alchemlyb/visualisation/mbar_matrix.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/alchemlyb/estimators/mbar_.py b/src/alchemlyb/estimators/mbar_.py index 6aa1821b..dd8127ae 100644 --- a/src/alchemlyb/estimators/mbar_.py +++ b/src/alchemlyb/estimators/mbar_.py @@ -45,7 +45,7 @@ class MBAR(BaseEstimator): states_ : list Lambda states for which free energy differences were obtained. - overlap_matrix: DataFrame + overlap_matrix: numpy.matrix The overlap matrix. """ diff --git a/src/alchemlyb/visualisation/mbar_matrix.py b/src/alchemlyb/visualisation/mbar_matrix.py index b5ee8c5d..682b8555 100644 --- a/src/alchemlyb/visualisation/mbar_matrix.py +++ b/src/alchemlyb/visualisation/mbar_matrix.py @@ -19,7 +19,7 @@ def plot_mbar_overlap_matrix(matrix, skip_lambda_index=[], ax=None): Parameters ---------- - matrix : DataFrame + matrix : numpy.matrix DataFrame of the overlap matrix obtained from :attr:`~alchemlyb.estimators.MBAR.overlap_matrix` skip_lambda_index : List @@ -38,7 +38,7 @@ def plot_mbar_overlap_matrix(matrix, skip_lambda_index=[], ax=None): ----- The code is taken and modified from : `Alchemical Analysis `_ - + ''' # Compute the size of the figure, if ax is not given. max_prob = matrix.max() From d4646f31fac00ed481c161e3abda6a4b156c992d Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 20:39:40 +0100 Subject: [PATCH 17/25] dhdl plot --- docs/images/dhdl_TI.png | Bin 0 -> 35345 bytes docs/visualisation.rst | 40 +++++ ...visualisation.plot_mbar_overlap_matrix.rst | 2 +- .../alchemlyb.visualisation.plot_ti_dhdl.rst | 22 +++ src/alchemlyb/estimators/mbar_.py | 3 - src/alchemlyb/estimators/ti_.py | 4 + src/alchemlyb/tests/test_visualisation.py | 22 ++- src/alchemlyb/visualisation/__init__.py | 3 +- src/alchemlyb/visualisation/ti_dhdl.py | 169 ++++++++++++++++++ 9 files changed, 258 insertions(+), 7 deletions(-) create mode 100644 docs/images/dhdl_TI.png create mode 100644 docs/visualisation/alchemlyb.visualisation.plot_ti_dhdl.rst create mode 100644 src/alchemlyb/visualisation/ti_dhdl.py diff --git a/docs/images/dhdl_TI.png b/docs/images/dhdl_TI.png new file mode 100644 index 0000000000000000000000000000000000000000..1e85be392ff773d7752b0c77dd94920d0525fb5b GIT binary patch literal 35345 zcmeFZ^;gwf+%3H6F6r(@1O(~cG=g+UN(e}Ir@*E`Kt$;f5Tv_9N>D)QPC>dm?%JN` z-tmrce|i6a%Q%k1V1Hx9C+D1NN4-!}z`>-zgg_uTN{UcT2n2}^0zuSAM+JZR=o#t^ z{wM4%r|+)iWbN){?rH^jZtm`E@8oW8`-;xf%GJ%*$?+-oQ*K@^IvaO)XE%{2PaOWw z1>8=qub(U|+@64kU^pupxIrLz=J0#5-J7agj%U^4Q2K8{$(xc$7 zV>4^%3qv-)?};T0LesMN1SsW5C|wvF46@FQLXa75=c=Mz5;=@qT|NYGX8w;NS{!>jx zLig_@OI|K182o48m3}b(zx#SCA@u()iZMr!`F9aX58~T@7jZ@WzX$!lH6kGR|7%B8 z40&t};^NeKY|7N(|E!ElKUhC0HWu64+nbJ&QFfc{pMhC-A%3F}KPRiSh#Y7CXKAYO zC@_K9dK{aXM_bbb(`_%e`Tt$X(gUs}{)sTtN^UcWy){cb-Jt$U5eiMEk7k(6ewpx6 zxS=2(Sgmd-xSx}a=lhzh#qasEI2KGuw6?1s_!1q+CyX($mb;PKsck&Dj92MPJ36wK}JQkbS-~6e=P-`fc zciK$dU?c=dNahHyVlc6iBp;AhSA!SZLe}~Ve?Xzy2)-J45Jg3$&f?U3A!rP2AOt+y z3I(MPW2fT|_1ibd-8BNMx^iYTta$!g&_h4H;Xn6=iL-G*1MkU;jVZ~5@Pz^lx3!^A zV$EOpUz<&*diDI;WgMl;zx@SnuB{(Do!l*B`B#%pgn~?rOvp>^-?h3z(4DQb{13Gr zj}~VyNpv&xLxybA=t;PB)S%F3Q7ABR4mJojs@WJrjy#Up){EtCloBlh0-b1S{$<_s z$HVF$61)W4V${6?q8F4}l9IC-J@^?Iz95H>BMii1vX(wvGV%mT9PZY1&N(5OoB#R) zK7=43c8m8AutcnQ?mO(N&PHV_i4v4y$d^{?+7OMSo70%%BTH`=gJr=D9cZ^NP0Gr` z#lyoZKyVgRKynqPL+?BpL0aoaqNo336AmjL2q{Gybyn-)$YhkF4ql0iAZ>YGG)2$+ z6vfN+pmSl6+M~{6jG-1~`)Vr8u&`Xy4+^z3;gTwY&0r z-1ov_Q_I5IC*>?IBGtlOT~hK9-krQxp?tOdI zpR3K`vOh0;dv~#<-{?%g81xX(dh%1$8e6R%WUExQ3bg8)kOT=t>!18T34At)JAZ3! zTSB2ax{Y(!al4cs>MbZF7Vmjy-!^n%#hG*xj;VJltuVRBRAZ2npKAqckIp%D#DQ1l z1>Rj>ZlsyM_-c5BD8b5#W*T_K?C8pma)%aaHWo&XEu~9x*o6in4Q(KLL6Hh@(E&@$ zm?Vd^`rLzxfECFF+We{5zN(#mCN+7>%%2&;VI^8c-TnP&rFs?Z6_*eQDJf|z`BPLv zRyD5Xy}md~?~PPFd~$LOR#sMXcL39P87s0;CShq?Mf?m}_zLo@WyU@TJZ54*xi z8p=+WN49jlz|f1yi~|J7oMf<>Y^uhh;qj>p?Im z>u>cr$(h>W^}lq!ytzs9DN9bS5oP=<^tYpoGq@ zQZ!~Yd8W3#xdtT2Y=e{Oet+vxAEl7rIagl-ZNl?U&&-e42lboX+3zpc@GM`ywirmH z=OM%X6RoQI76n}=s|zyUWHJ9-oPWI0SNB;EJrl!?rjV^TrzF-A%_RxWAAYGjr)`V# z_8G;IHD_&CRd#lENN{3eqAPX9@4uCngO4sRFPHZkmbg8a7!$-o=Im^YOQgAi0b*-S$K)DMB#9C=?OehWZm5v>?+TqXb>^0DOk3t zWeA1L&l~aE{z5J}XzO||G(9qDEcy$#O&LI725iY3#e<(qt zCk`1eA}kj9&oF|!2s{d`wjK?n+_zj+EQZ(zd}EFs!d zh3S)%>DE>;sDeVz%nZ5T#lD;;#)+Ju(;_yc*t8>X%)#WrG5B`Lge0sh7{m#kJ&eQx zFeM-2_$W0cga|yeyu9Fp4z%Nw!XI)_OZo7jt9#E|LI3s}H=`eo94-tDzrUEYiS)%& z$50A7=De9I7Ye$+4d`G0^XFAxEIEi?h}s;)*XMiikIjG?^~S$!^A(sZGYH$BC?!2B zgMjQGhd`Eg+cpK9OW8ib9?eUJbqJ)(awCK+KV0vJ0JBx7StmrQtAUQ_F_JFJy%)Kj~kfuzg&qN`0=9`Oj@H<_i?S$5{<{sWVmL@x2u!Qs=zzn`Boo7 z+Q3VCh}riP_UURXEQ#Ag+1=?XBw*gX;aH@HCnJ(Q6BER4n?D~@c%ndLBe1*1OSE{+ zyOEigm~ss3?4}!>Xde+0j(vVkX+Kjf_oYZp=Nk=vby>H}`mvSS06skLb|Fetq&*`g zNFo>F%+Sz#Nh;{(xC)%xY-J^Uel&46LY1WlOtaGKzyzC02>cYNwDR|FlyQqE=n!N=U#os*fEfvn+7%QP_4M>0*+S-Ny=yb^ zku5YdzkaOJEAVoX;!XoeB>n$N&3czEuN^bH@)%>jjC<TbLol?CN_j)sC*&W9?#-Tr=vJxvwbb8vHLP3T+X>h zi}3IXk#3KN7UrU0yTQrb7`%bJ1CIxKckgPB_U@gM-Vbd+qwJM5eIIQH8tZM7J&Qwu zVLbo9gAhtb)BFtDmY)Ws0SOVMtxSg{Bg zswruH(iE?~L@RR2?5o}TeRxsQ5%WxdiD9RX)hA(`z@7#RpHPd92pt`rt@-l$TA6yV z9z3N#9w9*zm){0bW;7`C7PcD-N(zUi*+dM!VU-V&PcHnUb_}eXw&kI}W&?d%{^v2G z0oVK>HsBKzhu`1btoua&uMk zLv;Wqu8BB0x2BS1fa|fF^C?B7sSdGT7hf(go>z;lKE*nvB;_V;(NTn~dBANHRaM{0 z=QkZ(g21u9u|aICuC6Y4PcQOOWlRF5hSqvcflrEddk`d5W8KFb5H;`x7(2llsaB!e zW%IR6>7?$ZOA|_yj$_LgnR(gonZ~7)4D*asX}`Evit;5_49^8>v7FFSWvu?BK($R+ zowA^{PW+p}+5~=5s>L9JY|YfLuCA_*Ho3hnGirL|wLjmbE_tWDdX&tn9uXO-FviT! zKi(bXpl_ZqFqYq0Y2G=I)F-*){;y_5Y*LL^tWa#^Fs|A@<+$~}R4B0itV|sPyZ2x^ zgN+mHBKk=U5B&a3DDY#p=O7qTkdb6WXJiDeL}rTUg7&tsITRdna3 z$Zk>2k^YeVD7*ZcBB?+((}4MEGg{mX4@PcfZN0iXQ-h9;&2u?)#!-;R5Z#9=LZ8X747~W!LIF};_W^5A|CQ(_Z}B1>)JINBO`cU4ePqgjaxvlV)fwfXFv~bXgJFIO&#e%fj-Tp>k)&>%02C8&^x|; z&uz@%fqeZW^}Em$CMIM>bu?UR#Pb8R+q)2w$FXe;0(KC-czNJnpWl_@uaSm zl0_8dCe3XdTvl&+9=>C~f_U>s6s$P!6j&t>8$r=(KKPDVzt!6%{lw_^@81OY_>ffP za8PsR0apNtN@eU0GaPY$e;*wON3JIVXRPM6y8HGxMzx9SxXy$&D)q!m4ke5MZ|}a7m)7)`inob9}}K7{LBe5dQf(l zRZLAwN(;fx3?QQ05OCjm!DUd5ED17dI?BL1W|hR`+;R;a=d*8|-Q2=J zcaSOS<$QbT!ev;C4lYJAthEVYR(Vgy#f2SQrD!|xg;=-3Bo6cqxdt`Xy|dG`wuH#i z!C<(pK2NVt|GF*RK}OpHZaBm@-5N&dA3yo3yqVJiU@?Tg%uOW3*tA{y5p{YLOy0YgC=X+s>Ur`~oZ>Dg; zPW6E>s>~%dTWcHE6NzW>`}6Y)Vmja;+#(7(yK0NZh5cGC89pkfrsd^Qb#_t<&3&O3 zzgf(-wL#GvB6$ya_byiD?h%0jMxAmTSqGWXkJjjP(buoWziTvRy1m@iw=6=I?o6bY zHRB=U0nzYLg0Sx?+xhNHSYx9Qzuov&$mZ$US#NhZ)>wOhUxttiqG`aH&R9TrudF$* z#@3%7pFlu|ws`Ixj>_T@pLr7Tkzr3advIoCW{!al15;vX^_=M-M6-`R{&P3t+hTq` z^5xCfyH_A);Et`bL_#eK*ZG9)uE(Kj(zxx57XsS>T6mF~xJuM9yW`7UVf4v^JWh*Z zFy>n}J^c^24Hx-beCgOuZD%A9vzcluxduT1*8l^aZVwPWJ~@eTTmt(E)qsLl1zt>L zNCZ#}4ti}?f;wyYd~eQXu~j(Zi%_U}{|5Dd;JW0c6bC6lLpL-U2hHudmj77KFFW^P zo85XZ7i9;#ZBBE_Ky|z1Y^mO4r)0flSrX`Kee3fJm7#@*f+AWiQ)C6K3)otVsvAs9 z%u44K86ja|^Uh!dViBjI_m6X5uz^M-I23WswZXNC3QznLzq+~_)Mtx%7Vsh&^nCEH z%3-z^A^24r-Dirx!rwcIiSJ$%y1j4#&@4xub8?S}NAUJU$gt`fOC%vHqSj}~%LWR) zNd|^@AJ7Rz+Hhq9X+v{^u;tVO{T30#~)G=x~c4Gv=!yezf8 zzHWB3+6xIzNg*mTtiv=nH{Tk~l?@+BPfLpgZX&~-M!F+jV%i}gCN5rO(T}_3-ki4> z*PW{u9wZTTFUD=s3JXKSTm6|SY6D_K%J(nNX_k^((OlD8K>%0<0U%v-g(b48s_#x# zEepB*epX0t{pKKqQ6W~eijcrSKO`3dGNxG7E=H`eAPEX;Ac@8V@wu!2!`%fW7zK+Y zCnpD9@3YtEn)x^^wNs;@pul^hs^@4{eB^Mi&d8wJU2LQL^vof&rpB$x%I9E@nRnul z6Ef{-OZm|Df$Rc-o-eKXnXvyB{n|Z)p(B|*^u1PN%*iSB@2a(9+}rt+*T{B=%xZG= zWxyoi?YrmRY&bkv?9VsJ4V~(yfaqADlnlxFhJ!m-_a<^F=z$t6O&@UHm_I>8(FqBW z^7&mL&P4x7(5{0h4L$FfnWJj9vD*c2ZUCFIGb!+syvd*Z!&5ZxoyTUAT+5vZQboxk z%ml*t1Pq?62lBgR`u6RW2HQfGYnXo0^08z9(@e+*Js&zAUd5=e<9uTX0urhjXvd1R zif~Q*k9u$7*|}KWKYkejE+jqkCVDrt0pEaUj`J{(OJwrLa?E=oQvEa^X1%OTQ)$$h zE%4Bjm;oMT!vhg?LZlQ#d|6M1l|1=1u8pCqtIvlS3bocYW9!_TpR;qTwSLY}*1^^ObCpsTO7pLwM6 zo-3rZbSEGoCN3@5C;$Al z;6>wLZZ<8z(FzL41ci;Nx`l!YeB5$O#5;TCwhN$^{n&Xfh3zhsedcq84wpw!v@08{ z%q|>rc`m;csk-gHJ*d0kV3o zUWLilbd^Pkn*y@56bN-k02*5qe8?YdOaAqo_rb#J*M44Zj3wyrVFkKiP;7_|PLzVe z^Hn|5gukkpM~jn=0?(uIlgPh{)L@>qmqZc==;lb9A(gC6$F7eJ;z0Q{}Y$r8S64B8hQ2xm95s-7^LSJcDz7 zeYB^iF-Xyl<8Y?}bPR9CA%=B`t7~!MtWYRR?bZOpZc8;D=XIIX^X!ixa8x%g47TK^ zuqQbUq=`f?0Gi8^lL&LIHF4Q&?L}EKYe-Kdf~jeatyYc_;hNg6j?qyV44p^qg=D}E z|LW_0i@<~{@8+MsH5C}DiveV9;Qs-mqXd5@?*UB0mATT98lTw1OW4%WRS~KK>@j#T zM`0w}Z*{`D;+1kF`-|_Xn*9uXqHWS6;!?vfFoVZn20V0-fuE?{#(a*hTWc0vG9RA_ zI0hGMDrKg??6p@O6=d&U6Nz#kbXK$l%2oPRek*|u?*RS*G^?p15eu)t6rS9-46a1r zz?WWC{w?fI?}E`&mn(4p>M2dUul}|o)Vc>_>rXGUJ8btk2xlHU6c>0+feVd;D=i0b zV0L#AXMG(W+wEU~fM)AsyVaP(C~>*(;Y&eU{{NFv?cQ0k4@jUircF)lAlTU5yaw0~#htHQjX%k-yKJ-`sLBo_{K|7)Pl zr$8_krbzwboe@&?dyd@3qT? zdo>nXInDB+%u@@%XZD%kwJlz5?o5TCUeM z%cw#SI&h>ub8CTx)2Rk-Zi(2z@QKjWn4#IUt7FXs2i0s-ZTV|M@by{T(4E9ZhKt~N zEx!`uqogx~ca6wnrt)islL96A_&T%u>kIgu#Ed_p0C}4FoG30F9b(jiI-HFsc2`^F zRp9XkUa1WboX6u6Kw*kM75mcs+Sl<(stv$?2b+wmQ1_IB)R-5|zDqr}6US*6M~nBZ z7w{Bt_G2gwiJbhGH4}r>agN&&8J$x^PneD1(Mw5p;SZd;K9Q4O@AayV)5QV<%SuypFcB1$om^V3b1Kdi$Fz9!A%F0k)$l=z0I;AMPxVJnUY0~4j4 z7fY0a5~s(Avvo)ydI<$)v@t?G74ebBk4wTqi>B{tiE)Y=U&K0B$!SZvR6)`SNM^3J z#l@tcdxXAt0moP}3d~qC1OS61rZPLoGw_iY8bzUxpxcXe5KxXdrHJ5;A%U6zURt8x z^H@8AKsRu_zela4hT48d;(hDV`I~TS6ybv~p^!@`c=xlPJu^s-9D-b;&Dg9?8Q^}6 zbhRb1IBaTCjVfeN0W!6Lb$aH2Xc9yXcr8yc`+LfPpYNYB66h$y{Lp{7{pNjQ zPUE;a{~p|3B80Hu3JfyIqp87*i5&x;8h3du4it;jI$via6sfCFB*DI`gBhQGjUp6O zHycZmW8)!Gq>d_vwj_JWvTG?*_t3-H_8r}6(RbN+pw@|{!5J^ZcpHiIa6?;ZN1T|< z2xiWaEc00^q$lG$YLn44XJi|(`_=`5U^b{@nGoO#)$1#sNU-YU=KfiP9NS8B4iPe(8ZKVLeEXLEs{oK?_^i|6)MO z4KWSEao@9>_bPa!w_I8Q%4ViG#q}+anQ*cxggkbLNN7x8{^F4>i5?(wqah-fSO{T* zCAZAR^iw!K4z7C>cOv&i=c$QFz=Zcp$MRb8D5^v-X@Ls6XL|@kLS+*C#$c)- zUKV#(_vd4yMmK4&T;X8W!u4RWe`X*^JWy>-k2O~7HaPaB)^nSr@G~T}eOWO^+H}nB z36;gl`Mm?b+s^6Oq=Zq33a+@gXm0SboZv9~h7K6NVjV0p0ZjddcsB)@DwI1pzuC(J zjy+4N-6K>T_Cpx3W|LbZPb9fi={+C}45furOLr$FT7A^hpm+oJO~>)Qk^YrU@zLrY z4!V(acyXRlYWbntSr0AigH*w5CJI!66p=^5AB>e+%Zqn?$OjhjkUE+HhGK@JiO%4oJV`X~r+J_JVukn7L@#qPW`hQi&8?ca_kH3N4XJ z6BCP-tru#_r23FOiw9DszPlv>z$Ec^X}CO?sx3sH99SOFBfxk5Jxelq3=SUNnui8I z_ui-XsJlz1{;n+@k}ryWE)4ETOJ~fhc|Q5hiK$RQ-Fqv0a%5Mdx02_cB=Tq`*FQi4 z_cK+BNQ5lN^H3%MX^4?gY+Co=9?%7tzf@A{98B(Gko)?~G4XO~aBqMAT^}o-Ksuc) z@V#j-E80Neb1c+vILrR{d`%$m4*cF!0&%Ab!dq3*evJvx1X*QeOn^d1K54&1eEyL? zNrZ%8h9W%gba3zC7Ri3r5^MnyItGBA#JOUV=LWgqaBLYy-aA`R#@b$ilP`rSbdzZsv2u07FAC$SV8wpbI)@1Pvl#!IF2#8bnuN zEf_{kt|(nyU7&+N`zYvS;(h~utFjs*0e$FAd%r56H~Ok7GpSlzRM*6HT)JjGzlzFp`5A2Z!(52iJmEkf z4587br*T$S^PS73E(jGU<+c24VuW#hLa{1e%MaH*d08pvi1m_YI*O_`zlK(y)DSF_ zK)yDwRy46i3#Is?+ZE2Pc4$WD!OCvA#!o6(mFCW`TOXPZ@#;-xF~e06)L)Ohx&Uao z`NrQ=7N;5pnvx@Q;D#KxjT@}**L){E>K(f3H+Qpt3&sMH%@n*5zBcA7wZ zQCJD#B>;GB1P?hp@1hN*1?HV;8Uw| z%a<-|XXDyyP_M#cN9on=?CzZ6WP?Rqf(> zGgH%-$HdyNA8f<#<>kxDM-O1et}`@??eh|$piQy7*1Upv=(EqAFCc=cl&iOE1@Wtu zaHRiWvCZY~+N1SmzeQ0=N&4>ktjhPy4k!^+sN-M*?#Nw2**!y8uoxto0HT2*q<6`2 zD1%ogeJ1Y!ld7+ylq$gDW)=i$A6M{+)xiM7?nl+q& z0*eCjlxW#<*uasHF!=>VR5apjEr@4#$5v~%N!@MAFx3Jp8gB~dr3el6~@av17agObi~ad#aQ@j31(^m3q8hBn>x-rfNJ@-I!E`5K`{KYnX;h! zE3}1X57`$lfKaM4nwYa7-tvq^Ju4RA$+G!lQl%4l3JEkqKF0;)Rb{hf1~vGkq-a1B z09c`XpOZiOb#|%aRrcTeX@jt(@ejXqR;27)Vi2-O;lF*Awzai|S(l!X5eq<@PM}f& zSUW5(t{0AvC$+x)c=(PxtyDlgyN+lw;sMW+Xi)*qrl^kdOe zULMs1sBYdyM)rV@9vjwnF1Gpo_@%3qeHTi%+!Qe)&8&yO5F7{jjC6_L=?@WnK7NS; ziFpveK!)_kECtw3F>{gCK(t!ad%t-SCb_coLUAM;_hg4~SB9ibSM$qh&6OM2FOTYI zliaxm85|BSU!k$of70KrcHWhVsckzf3BRZluUdfGn>M$7`)#$&m+SJYb7`-OGeu&{ zPu|Xs0~L=ZHYbO!vAH=AM@Zi1K`ADlw?$rk1r&zCoCM)3FKS=%-iLIr zUcD-ae|>hozYqiHc1G^qLF!0 z2D3uE5vP&O5Q38+DLVFIw84_W-jJdFjuyHGh`C1({_kg~L23^*(Q!Iu4i|J@<3UXl zt%z-j*=Y&J9lPD#Rai$I;SXW6HPP7zcs7)fFM=iU`%dmpnz{A9mrhy>&joI`6ydKB zYq;N??*G2Dn58W~d6D@@J*{L2S{77Rr5_Mn-ToS^9F4fIREhQ~$mNkOEuz-LA2d$e zH-of>`VimR>rRy$kJY`=1DsjPM5!KWiv*wtfJAB+Xwd+|o(L@gN(O-D8E5m`J_3MM zI1nsUL41OThgDAxmS~p(Zyp16Mh~XAZte>@x*S(orq8PBp`)WOhBHMQw(18429z^P z^KS2J<`02u!~#KEE>Lk;JNfm;Q*QyFnmF*1c4l^7CyQiGSbYO4dK#`=ehloaGW zmqF-emS4G%gAfO97lW&rtu2tGJ2^Q;B~PQanIsgdq%m-EV*U6em?jcGK`HFY1fQ%s zz|*@E#Q!gOz((LMZP0Cqd>lnCkeZCz)xi=H`Wu!4U#?*56awiLU_-=ic5CpFyLR{L zY{n2NpV@Z;(HjH;%BSvnhq>W<+2_*I2&cZ{fGdhw*azxCAcRP-+j}@3K11e1$ek6h zXm}co4`HE*v&1^R^4)!SBbhO-hFDMtT#OHwM2;@HDpFUB(e_|Rh;-S31$eIrpXSFi~DfJ`TXsxZ*@FzjM zI=hK(V1c|Yt585}Pyl`2s0}}!|1AJ-;cgjABTnJee#HWJ`ZNJdxD=U6$%LcIOTW7wiANeuvZY~7hv{1yMg3ag$@iZ93g`TA) zS^!RW0jffy&DW!*{?3%0y-%1l?^^Q2WE{fxm_qJRk!M!-lY@uQ+g5)Ye||Ed06HV_ z;SLe~7ykGk-E9y@q}1ozT`4Ks=Yr+dNX)m=Wa*P%?l%3DB7>R#QR#{>8>i~^VmP%d4mJQp7S zqJxEjfuX3P((|O_HWW_Mf{1BOpcy@s_GAsnC;OT`cJ^;{TE+J7dL_8Iae-sWd3y^5 zUjJpqrWQF3Nf0E%wiz!bkhtDrgut0)@Zcj5W%X+zSNAsuaJNTBK`|(3%MQH8h#`60 zz1kZ?Ku8$&_APSj@pn#iOiaF;6O&3a1bEO!MnjY3PRo}LL*1HdcyfNQSSj{_T;^ z!-xz0ip)f0kWOcp4jRE6Jqf1HpsCYSk923Prt|`c21n{8&=~RnvfA6v5B^%vBy7#r zmCBI>R(eo{%B=?jKM4n60t@gwGKz|@ahI|z3;yc_M9JaRX0~kZ3!+a2#xz-3S;Yo5EI>!9RlbIVgml1(3kZGV zTyrBGkE2z&GSiO5ytU5$Sn@EJwZ2^SYzcVw2G)DDU0U%92<5DuTHvBHpj6{~`qbs> z*b2x-(5|koVA0j{Y3dNUlL#6-A1c?CETvHc$S{w1*KDzo_95=MykaK4JMA1bUlTT( zE={F3ZUx4RzUa)$`oyqo6D(Sd#z|D(bZl z7d+_$jh^FTtL@Q6+hfqvu1i5Rz7=ca!|6l&`NpwVgl>CrU1n?)4=NXaKzZk~K0pB2 zTOwl$0lT+An#4#?kKpRltJ^;~Jp2u0Qvx#mc9*%wFY8;KeEw<`snLsxwaV=2_RcrC zHwAZ=&^NJ~kEfLNJ;j0`r{;AHXRqX3BNH(|7vGJvKBLMo}P zA1}6ggRmBi$Tow7mYp8wvB_kcR)sDR{i-Np6cwkW$VS6U=T8me#AB$mHcNLisXnfG zu`qYqgD5U07Ogu6o|^qH!pK$*SOZB8$&K9U?bV5*q9V|*2zGaOgJ!7>lul8-F{Eqo z6epk0PVr=KWd#B7!BM{rijLO$^TBSSV`0gH0*d5`Nh}Bg@PTS7)EH94{b~1+v#X`G+px=@6T|0 zn6TR)?d>aVI4cO&=xF%qh_=nu@|W$&uaqp$KfVPX{|zY9Ky3QJYk>;)q?(2{>mL*# zF9Vxd0U?0wsWo<9em<||0Df}ASi%vjzarLYx3NE(@dQoKN@}Iqk-CFxy z<0{g>Q@=He;xHnj-%r$IRR-M$3VA%M58%8anH({}J1}$7eq+==@tI?xezPOea${xs zy~Wn)t0F#i=4ggcfellVYw3f!?4Wz{z`#JTls!p|3P42rBp=koTjM3WRNp`?0Yn2g zqd;Mlw*&r+3K2UQe)>TWNXAOa=0N&}KtMiv?0Ewg>G3m-$w|0JhN0rXV?Izqg#dj7 zxCj-LV<{pZ7nGHiWtQ3lTu%Nxzr4E21tr{c-J7WZmHGFdfeY%l`wPQcXRzJ0z|1aB zf2;qHAa5B3a%-2rznEkEk27JbkdNb$+<}aKmy*v z9wN~2X_q|zBJvXferByD|6M6h>&oG3^4Op`DD_u1c2S&2_x-GJ;Fusmk}<`kGZ}W? z2FptF+4pg)JOH@H*kdZ=fua)hVY`}dUn4hbA9S1GhXTs3SfO~Fcgys@7R@ovMr?Ke zUF&+{DXO?-=r9u{{yW?M4$8Gm)cNWm+j7-ViW?LGfTCE(j)V>_T!;5=kAYl}AWmM2 z<&ZNWAH1>~978dz!+6iclDymj&ImD7zmcWH4u_~ocxr=i3H6qjqF+1rnhgE#_Bu!d z|L5~V&m;Nl<^Un`g%C?p%lwae0worsLH`XYk~fclHD44{g9w3O9`*|} zeVsM9V14Ja(2oF_(7Rz*6EI;{)#ss~a?0zu<4r^{3G&8W24aE5XuA?`2IWe#vk47%!vhH7B*xZRd{&yyhYWlZFRv$R-mH0l5#e;~iO4XP}EJJ*GXPeW5h$BP}6JLYq zXf5yT(tMvILo8E!mAB6Bt=U#$IoWgv5>N<(SI!7?X7HAE@eI!k5c+C7@u1?%1rGtE zq)QtFfC4km4St!Z7vIsP*1>9=U`OxZZ>2)+7QbWnM}6&Ak7TQlsJcT*QX$m4)dp$quiJ43UPUroN4ek9S9rs}ie(K8&#B+z13A%q` z6~G3(O?oc>YArlYv!9M|@aITG>915281d5`{IVrodlX2mjV_vL(7k3vX{;F#xv%=d&(L@4!+O^1MNC6*>|_M)3A@ zo^qA6-ojt2QiwO<(c%EN<+Y@^GMm_3BCaE{a?ORswqn2aecKIA3xSOaH^55lF7)i_ z`TkHUvr-kFzznyzhD{2wesSqtIHB~PGg^T+td_f=Kne;LIW8%0B81{|NY z>eSFZHt}-mxwO(S05pdvu!~S-ZEYOkW!4dSbT_XfkOnG{@3 z_4L*K#(%W?dLQnq(2FCjo0Y$_>az4tcbsZ1xf%XZGnl@RnA;uA77>eX(XhimlmRYV z1HwSqNZ#;ZY&|vf8)s}Rx0sX?B~{rrN(la#LJF|Z?`eQ^82bxdrU_b;>Q9WVXPtQ) zZOVnvwTQaWM)Sui&s4yb+Kykl4RhbCKh9|N(cEU!MgYpNHDq`0t~Y1H638)dxbHnU zwx#i0UOJf0gXl{zwlQj}aqG(7>SGTgU&z1s41M6`aiq&wR|B)f8~{NW&_bb|-D9%` zmgu_DAv8PBs5;BV`(2J+aq!i+nvJHlHNT}pzx0^#dk=g42sm-0_tC4WIXDmq3J;|e z@KPAZ5MaOWRz$5a{=jCSOFBPk=z zv*;AlyK{5UkpN^^#Ku!L{Et3^IGXo`Y|I8ZQ0OuYe&;Fb8tyrWbCkn|!BoP&Gy$lp zq1+viJai#Zb^>y8bUZJ;0n2m2YnA*lI=Wyde_*wKqddj*+i2tcD=9rZ!8f7$H6Uy5 z2Rq7wN2q=}dV1`R;;hNebJX^Fv_Vf)NVIf#rl3IL{?fu9oUenI#2tZmhoFm*$z{7^ z{y8)xYhq&ZPh79^jR_3T)ps-3>u~?G*N}%j;OuMd9Kz>{H?*)MxFv`KsphyB>$6Ku z3tba8((5Q+jr@y0scBi=iPD*S&6odFxbA;7?PqH-fPAiBfi~bRD6HOJeHct;1!r9H zw;oYY*yYy3cbdldi&seLFP&|&Q^?HK4Ei{-E0PkW`D91oycjhkCt&XSkA0s#fMO%r zwLYi#41=Pc0JPI^?qq=k0I>kp2i2Vg=*65}IeduYz~`iXM!voxG4MnBoU32l2t{t7*U$tl!&kJC|+u zN5@$9ZfCG56F(dvY~()qP+VO6$jJ#1aX_Q$ywZ&V+VtbT7&SGuLJ}a6f?86W*qkD^ zPoIRt1O949DJ6<+D3hAZhSelu!JExKArBD@Q7`I-!G)2#+gl*U*DL8~VPcYENtQ0% z(1(K2zSQlnz%>l7?~X5^=7D9zIy$&wx-RFruVCsw`&Q;A6@e3gd~f(<_~RWKw{g4m z&(BPCZkuXwazV&_3m?E4k)Tak0o_V(Fapx@#s;&873i#{Q$r}6ChLu-3nZ2 z4`zv|J{kDPr|f1YP@G zpnSsw`JvkDK?$4oiv84g?Zen ziFm21y;ro)`an`09tjLLaSCo};j|98e0bmv6O&~QpD7=H1wWhgcyWx}@)9u{;mDI)Pbn3Kq|+{r5c?iQ!(%N&VI zdoT$)41f3lPhP$_<&qyxq?kwwFy4tE>Vz6ROtjwHsk>!mfrlybNt8=$QU9jIO5jQn ziS9XeIZOI?5mq8^zQMm<5u<$Xi;Dw1)s#!I?epEmA2l!!qz$)b$C_WuOET;!r zm{>yRNVIH8=~qq<#PW_Cz8c z_HWgwDf1jN(X|BN8Ip-zvW|9Eu2xprI(l4#3a+#$dj={oZ6y{?DMZ9pUK=R|Y>H6= zVPT)k*K-7A#vQaflLrdhb)eGR;(COe-8-)So}O(Kjhe*(oOn~}5{nni=~l#Cq@=9w`_t%|qR$W^R3iQ)ZzlKpn{jSj zE(d&=Oq-%;Ny~c@u%=v!;g#@AtJ#Ga85ATrst9y~x_~$EjUHrZ*&L54jqr7u;5n}| zQB-Qu;n#jMcMip!fZNumBAqt~5Y#b9y4hqv)?ZUlPZ^UewVy2J8a6eF+}7-&PnON; zt-ko$32}1Yo<<+R-hzqL^*p9dg&UYrt=(fIBJP2~kr254Ab#P*{rqkAtU7#YNS;XT zZ|s)I6yf5GmB`9M`f?55Owo4Ep4~Wu!5cOCVx5qdg=bfw0bvekDh&h61D$o68AK&P z?S$fdnd*=X@KbPf%+Nyq0;>G|ahU;FRl(iJn3!a6o(cRY7wVR*6rIZgeeP|i>7O4; z>{<*zCDp*YQFL{J=Q-&*DdD4ijp`06Qm{N@NOHdNtLt%q6sz>})ffTcC%&ct{v=g8 zMZ#|jqmKcN3BGGWDETy)iP==J*LI)q(J(h$Pf(2$XeTMPtgQdOQ<$#)2NQ~7{vXEY zE+1ckt`xYNY9PE#WN;h(QYS|^wxI0@Xt$Skg%sd?4cJHG4+RijbtlG_e-)yq&4gib zsvTfO0imn`4miTUFaivjh5i52_f;GG#C#1d{3qdJ^HibNaz^12HG79xpgvOnSqKLn zcb>-_#&(YFD1cLF4d7b|>h6RddL~^Wh~S?^A=`!uMbF1JBEAnW`}dVNmJEh?3lNFty z-bLpRUz4D^qVQkf`2;$d*`>7O90CAD%z>-G`}fBWOS;5$cZ20H6ALAs`cBg&8QDu- zQIb^u`(~b(s|C(7IV|`^jv4QOd^rir6|RW@ohm{NDIM)|gV>P{Lyht_8RRY_zi9Y( z`WVO;C@|p`QoRGJ>(yJl8?#A?vO6@jwVEe-5Xh_;QDoek&*pdA?2y=A`R|2RLV=fl zTn2|MB3^dL8pt1?MfD{ixh48uvHsUM@L>UMSovEkRIuyJ-g|FBFf2sNgLW;Di>7N5q;eA_4=?G%L_q)badrAA=EKdJ;r|}MS3|5iGT(;o42<&LF8Lki!Z0P&E0 z&#kh}=2-~v3pPuVmz}A@10cIehqw^Zn3^W!&5he_=1DHKrUx!rgdzrjFK6H``0$hh zZVm=)8_K{U5>rPIziES<;BxN*??ib)U$?f03mF(5U~m30^1o_(51_1;u50ii34#cS zf=JF81SIDuNf4Bz1Q7%kkQ^iNkXCbxPgZ6s?vYLRwiF^Y*=?a3y?f))cDHEqYfl`kuxI(Qj9bq@8|eq%io5s{H33p zUpO}) zKwd@}+)OOrqtG@F7L_P{RTbA8UuWsyDQrHig-L=x2`5cWkqxp)mFxUB&Ja@ieF^(%*M|T=ag&aaUV4 zX_rs?!1I`j7Wr)P3wGm#02Er0y9Lr;ScJ{rZ(-BoNNqpk{NY$ZgCpk#%JYme6g?UlsoyQ7D)pn3cpJe&Y^zpxyA1Dc;1`J>(3qI@j%9h4i6Wks)o|BjL4cWjB=!$>S z(?hX#qD)N(5lYvsrI!GR=QZpDy#%R3&pC)7h6)?AorrM>=$KdT39-V&=@G+z5TlCq zuIfPYyQZx^TI*4K5$Z`_s(>63>>AgmXC+GoX(X#r_FnoJ(kR%_1k4L@6NGH(ko2J8 zd#uEabJQ}NJym9r?3%RI;Hz_(e13BL~Kp0qQe1VQDjqzZ*!Mt19{M8t&!CB-<(E6Q&Zed>4_G8gp%~l`}fbbzQEL(3WJAF#>{;@vNE#o z6S|CNIjq)G<_7W9o{hN1q2VpZP5*DvZjVx;Vqgh$r)qGDt-0WF9GhSxkPS9*#{;SA zi;w#l${!m=_Vf%7lsUjto<#`AhWz3A6KbXjhAEA~1St9bJDQrx9>1qs zqLGSrff<68%{OC`q4)Zf#y-CVQ<8^km5Sr(_GrZxjzjF$^`MiacUpt#QQn{87+yvZImZm;8~K`M+)|0} zw053!4e$I7)Y@^~qJ|!5zih(7-hD8Q7cn(hYe0kZqReEXsNb+2>R7E}Ol)kQ^=EN$ zaBAL*y9JNLzIk)c!RQI3jN{XH*Sx+W11@^ddjvw>WR`o=qT!VgJKCU^zvFd7L%-@Q zL@n@-NNn3uhk>fnoBrjj?S}R_R87-;-+=+eo@`>{F0a>6WLn+W;A<%ZK4rL4G}8ps z>w(y5LermpA0b`!WY&8jU{bg$>=ukVwCYLN=Kcjeww3`nQ7WrLtRa_a{DClUIid*L zTA^uymJ^hPewRDq0Z8`3F!!d647}i=3jvS9v@Kdln4IqsOP>_D24DV zaSHf~iyi9A*hQLQ90nwKk+=^Df;@N^c`KUK^*{eJkvx9tKmX!x*SD3bkewZuitSOi z=C(FDsI*T`P32jSNua#-xtff;4zbGg!$BTs0$SFqHy#rr#!4@vk<_(EEKv*u1ubbV z?9Hu%coI2U_dkYwiH0L&^iZ-Pt@Vo0rtzQI>d6$=koL$7J7@!2e6UN-vCrL>N;PM|3FDzre*Sw5)gra7B=t_}_fQ|`9|kMj{^@Ir#jy_`F&#`U`OsB5+0ddgCpt%OIy)WXx|AA0I`xa& z)0YK3`+jd@^6u1G`0x;d^u#M2C%UWPTZFXvs5}!mydRp-<&6gljxVLw`S_;u?36ue z?&L0ADcxq~Y#=qt^=3i0q7JUFE1Mdhtm{pID7jez30<|Gc3@DR->PYYrJ4fC9~^C)j8CT) z%AK5`Ay9FNdJh3v#llYSM9!L1+rz~MiA_0+)HHLGgW)lgwfr!lNSkxaYKAXZ&>@9t z*e)>g)ogiWfn>W9Sl_VgX|933YZD~>lZ8k~#Fwu~p!m@|X+BiYa0OpOPGhaB$lip6 z8kvJ_7)F#0;}>kG2=r`C9T2#gb)jRT6siv{zxVn}1kbwMWwztF_4w@1Vamq8*eIw; zhxcfFc_V!TiQvs&Ssy{p&USEB8g+E946ioon{t^Eq!7O&BVmM>e6_Rv?`0~)a^dAD z>j1HbpM++sok5%?`EfYjH)f>Qo-hGBnTJ9|hiH&^LHrOSaxNhXibhCq&Z>5TaP-rO z)7!59QOY#YEqWutUxpt0j!y}V*X+hgrj?9+9$>%MH_1Q+dKz-HwM|ClurU^O4Lb1ba-p_nzb4QSL4yq&H8~-Vk>k zr_|gAm`FO_lwQVuuxKB6I*;&jy=o=M}#C;bvZd|bDIw2P^cr}eAA>z8HKz^rgUO- zHY_udJv?8OLG^Oi4k5~-Gq6j6)Lr12a!tOGfXQ!J2BKQ$n+2SQ85nlcZEq>Zc9%uI?j-EeM5rHrC>&&*WPWtr>h8 zjYFzWi-WFHw-=`$356GoDAOxnY^}W_Qd1Y_eid+PAM%_J-0rO1>%D_iu7vCb$d27M zRc<_$xHRUzPYSqwYLFuGCI|o{lDO+UI1v<}Y-wIyrg`~twkzHIzWZJa4iwqYP(fiY z8@_NG){QP3P5K*0(DJd9u2de=U%T;~C|bfWw8x?(El3$@$qq^SVGmGgLQf^cBW$8o zLYthO^@$fDVjLyCI_G53bb*NWdxn>?vu*^(CIA(JibhpbW{~3^ayqfkgS_Nep2@ zHn+^)wnyP~`PInJ(>Y$qC(P0fDB%F-19b=Mm(1lQ4LKuCY2ha?)}^DTUO3}St>~By zdrg4BaVWl`iv01k@AWqms9{8TK>i=uVoE6AWEi1Cdj~4UInbRKmY+{@;C#}PC3?-Z zD=%)cU2r<3nA^+Li2;U>jBIeA{17KU-#+c}3Djl>P}@(TClpoKnB%y*;ol$HOK~w9 zG`%U>|LQtoxHRNJfraIo@fx?#QiJrJD8Ecxco-q!xK-v!;pY{SS2z*GiUeo{ZH2{( zpz8s~80jZ#p69x3h=5#hsG_4TlkC1gP0!A$_W~Bay`!qMyu0dRll~9Q4nqSXzyicc{I~-E9xJ7D3rK{`U7x<*B><+!d za?e%#Ns9;Dqr|$EX54)_V-cy-5e)3&_l>h5gKbe&I518EXns z08;80Yp2Fd=;c;^s-5Wl*F|+JoSN5Cq}WnCbjfShF(^4Nj=sv!af&cYEsV&3y>P%} zs4!giZJ%xY!NhH(jFZ*BRR6}&&6Z&e{nab$@xa}zwfckJpO5v|6yun7mZFJ11)pGf z!Vr3wwdHx4K0=#x0R-}kekJ6rshS@gvvf^XsDU(FkydCC}1cF+PD+oEkfO3G8a~@CenUoW%eQbQM1&EO&Ss>SVQ8-`kj} zJd2CFUR%)V@E!hD09d`xTFl!@H2e(Ki|1B_2rqVQh4S20OgMauZu_kbVpO3FWO z!eh56x3$8E3Fl|&w`8?t8?m9Z^Ea;bD+eOsyB+2l+8Uw)11mo`bAo^v@N~Z~Kso1? z2(HqHI;`gIB!Rd^A!*cT!gT_^^va+nswell%zkikb5pbHtkiZ^KdYUj9=ma06zSxj zOS9C0=(0ZWoJ=J&q@bCDbN=PsQk34&FSkdd*4Kkd9YrzRw%;`*r;)Es9!to|HQda& z&zzi3rc!9JZF}&qmE>S;Tg}vTT8e4P^b&M1{}zl4;e-(T2C&P3?gHZ`$ENRH_?QZK zjZ|J=g!i_d1$IP)HGLWMn#=;C*@Vrpoa>q}`@$-Meb5Pdw^PBlw zwxlFD78G&8?KLg@0S_3Y;3l9SJ%uH8_X?L38!|!&2>i{7I`^M2FK>R{WCW$|R;$_4 zuorGJPHDNw$NpbdKhw>J?MRVFmh1ZgoSoV?4;aiD0PS1KchG}P7D7jC2pus|pXIS0 zwwHTcaSa77-S=DPk(Ry~lu$nxnezMc=hor}=iK7Ct;c*B2JhO5;*gGATJu>eeP!GD zH}K)38S`2w%Z_r5YTl%M8DEKNX>SeG)H^amUhl}+bb8i#5jmpOepkB9A;4XhMH3=h zO>LYs->`Bb3;qT;+14-cWi5O1xBH)`S;vJM%$5KwS7O2oWg6bS6!tg)i*V} z9g#CL>nPr@et-4>Kz%OxH`x<~--aGn?lE5Dr+#T&b%>XsCv2wb*$*}$Ka=J0Ko3#( z>sLX838r-EN3Lm81@ljHistjfxYnPygc8U9f+l2(ZGMAM@>aGthA@Tp#9;O<8v&u6&E;&40)6>*75u=uqr}>L za{MD->2p3xjVK=%;{lTyiRz>z?1##F`GVJfViS*&`Ssu2Zm4S{UmYP}U)$ATzrH_Q zz4|Yj!1#M1DckzpA6AnMFHzpW_l|j%$2e5q4>n<6I1)Q2GjL0959HvaZTrLE(skC?mipJrwRM0Nj1VO$P5e%Lo%+M8&7FMwOg$0m>*-Ca=e;7q!=YP*)XI>k2Z z1hOPJRG+bxuoWp4w_B1Ey}1X~K3-n?@x`f&2}?(r?Co#*a(?dq&x*D|&PF9T?#g?G zM?~Dz)+PZL8ErE(-skJn$rA?18w4Dnflb~D;ay<}v1f3Iu+3jRa-qIUakr6RqklZ` zd6vzer7NiKR3zC+%x`Yw_`?T@HvW&R6PVx?^+FnnypWGIAtk1zq6 zTU&w2f9GWrLQelLh{8|h5+(T3VoqJ+Q2PGxyF0e91KMc5KJZB4K+FJ@&6NGAbq0&M z(~Mq(qdO;};1K+Yf2u|)%n_to`vwjwa`WoAz2U15gNZi20n*%7m#i#)XcM=3gq zRym5t+!lp0?_ip0@?Ttg+|RVLw-GnK^ra5C96S$Jn_P~8@sM*%oa5L-^iPM8uc3!z zO}@N!%JTw|7!Q$0zeuM&zx(ynD(z3D%5Sz$>lIab<#W{0xbHv7z z+T$M);@MX)jhsIF(<^@nxXrnLckr)}U}!c2Y3g3)A6(tLWgb$vA-*uywN;cx2vmFB z*ZvtlSy8Bam z#yj>A$!1}M6>F~6mcwGipiy#}%3fjHxcg%U*sJn@E+CeHreIJNxXpt9smvg*l(gntFF6ZONut{-f zP|U4|T0FqUawQ#!jdwxAU(n%pO3aPG{^NegruI$^=(&RI8xB<9^!1!jeyvIc>~YXg z2IUU5wZ7QEz4ZYh5YWjV?5s{OB@rOO=vHo(mj}zEKL$#7arJifAmTeaA-`yjMRICT z!0`fQ>P)^C`4TMb{fB^V(6Wn`4CEyw2%}J8-z@xx5bQKx0BQkFWn{*{1ArPm;G5bX zH0`DSa$a0G`eyOm*mkvd_Sko*w`^|U2MIPoQE(#^R>miCxXP`5FIMun_n2q&KNh-I z5-MMmfw)>3HztW9Ie%$p*=KavlLm`*k99NY;{$OpGcZ!skk+w+7=wN5t`z<4UC#T1 zgY)nEtPxie0LXcEn^euF=v2>0BNBcYeQCRcEQ5SD$Abkare?$q;3f|l@R-u{TdsV? zkpk?5HN^KjzJ+n^Z(Ef2qG0z2ZPEH(f0elTLkvE?xvJfL&g4eg_GW__UOt}Q`CwrWCbF@^ zuRl0unf^1*Cr%pmrt6!#r+K^H`DpywH>&Hc_K z5^=WlgPO=4`IA)RDR?B5dphQ=EJ$Rg+$GbgNj~PCXVc}k@cV$gr#JGkpZ!*?A9CE! za~+X5DA>FXh5O)O;sn)z)$5#G`l{gKC|^nOw1{u0C^N8wlKkUTJ^uJi8qXtymqYiK zi27=TUk~QnBqSX52W=|=Q2-uj*QH`LSGs$M{mfN&`?vUu3d-J=gRWOC`u^4rL81N7 z^;*X}K{Z{13#$t0;rLx?CiI0*^+?&Y=IaaLg%6^TPxJyx8N|__IaK2O6+0Cn9K6DF zAI*;2-294&^NovXFq8WPiQfTb>nc>bvCybhAcuL@S`|=N04>^m%Y5F9!7jU#b3C__JL0TE6+(MZBG)Y zs|Xy6H5Z=YCYd|lE=c&V@dVg>j}w82lyP=}iuULCUgA(dq~@{a^=E&;I-_qI^!MRu zkvI()xRHp(ltS9f^eoM9K7A;Sozkp{^fFHt+}$t~X1Oij6E@IM!D5Mu`3!hJVsR3x z+%@kA;YZ@J^!90TWSRzSqo+GPkaKjvkUIsdzQohu2=wD2KjAq0 ziTkXfEQ@~=d-X5Ya&+UoF3~bH&UAmDVgHTAQRY(sV0@K`JeW{^z?*Tq`{@Pd#olkd#U0s?&)U) zXe4N;5Yrj4+Yl!ste>(yB})GZVJv}rn7z~g8>_hDQ>38Yc`cQEOm2r7hyM~ zqgT&9OOzhwD}6(eGUYCuQecnmO{05kSIaWFGLjh9gj!5(5RVbD(n0 zB+MAtP5nTuKd*{7u+Y)dpIv-0QMopB25en54_=n zd;rRZ8in@ozztJa!T_5B%mv9G`Nv2RZ%!fGNlZX3%{t$Qfk=)Xq9~Bf$s!>&)Jaf^ z>gqS4@B`MC*od!+kli#6BF%Dfb)5uVVB8KFoTgEgPY4H4Z~+4#*}yV@GR=92`#>`h z_5e-uB7{KFP)!QXxhXF{0ad^=}+`+DrC>lmatf;6)XU0vWD%#n)c=W90zq*3o zd5CNkg0NM{Z%84Q9{kuTK(5{?-MPzp6KuWX#?x?Z=cpoq=PXU=$IH#EmdLUdKRs>g zFyD9RPrXFpl*=OE^^rh+WBvk*7F2xP-DOzfl2dn%sW zS8Vp(Z|>|2)Xc<0ArU=DjoEB>uAoj_Nirv?Ns!f!y4|14C15nb@#GxPG_wYy{VaWD z32NDLVFaMwMgmr-@e&>)IY!mtV2+xa`uni)12~VuR3YmYY%PP~vWSU-P6_tD`jaz3 zxvQ#6^QKfJli~u7sR>0nc4`H3V*sCd$2iOy&Ux0VAHFQXaGMzc!54#XBp{WRgC#n# zxESI*<`Ji|io}b9g%_>yK~M17F3ptepl62smdss6LcD6}=;`iXzYyIVVkh(M8+c*B zOjy?Z0~6EjABJ#4j=@Az$Ef7DANEqQVnvwXx1HH7+?;1~O3%QzxH>?HwTBLzJci{! zH$Rnbu0gJ|z3K~){J^eb(kXmKP~u{*Uoy}mLrPSN!_M}YU!RMRN_u}m?)^f8GMc6u zR}5O57xKf5!)H5T04R(CeDI@kJ4@!Sb5v&B92=-2v4~hY<j1qOR5>=x+u@Bv zUw};x=n`ii)gZxH2e`mR3JNjfuCb2-j9ajT!EfoS;(b%Ev;7J#1o$ab%Cr9O`PYVv zcu@CAbI$BO`<;)7r{nm9X@art=HVjZ0aj2VyM=fK!cbCz0tC@0nHk{&BrkLj{P|Q@ z7DdTI5fF1f(xT5;_gt5^Ax~wSlV|=}-DpW};`*1uiW>!qM(4=h1l~ZDYJYq1`4tR7-gnc*+281K7*4CA?6L-r+$eqC<^fFUBTeeOqg5v;@K}#VT4P; zd|aG+By~3CJa1vhW9;mE7Tp`Krq1G%n{fCN&67P@@BMlq*{43Ync?gKiI9iS;DJN_ zmdFKQ-;KW6*%O_e`Fy<6X5GfRYIe+GhKshu4D@HAR~#7HBc%}8IE}OoI8YsOdnupB z49T@Lcb$Dw7B+R9ZJgZoPkvYd=y2DEMe^OVP-f|sYUD>H8)DNLbFC*=W1)l?w_u5jxG65id|1ZK1b9WWxU$N;XKBfd3a`h!6ZW?K(@9IVnO+JVyGP^`{{Q z8eHZS@L+o!`{K#$WezvVWHz5opB)O2*)U&nNt>(M07OBQQgoZm?rat)Ie(ahtai_j6 zv$^Z^47a0!_x0netjOr3J4;TWPpT>Pxs*nCv+lH$K9}+9_d1`Ry-Ezq^5JK%8HPrS z)?F=BTsITvx+ozH!kKJasp&vqQKf)mHS5wCbrngks>H~pNaI|N7tQ4EcP12F%7v4kzdgT zAvf;52>rj5#k z!&mrQv)U{?H?&1axJUdv_SHwXC02fAaYD-57OV5a%fxbc-AZ)9?beT3&PT(hu1jc{ zesFjqM%Ktd3^rcpP+<0P&hudHCgL;=M)rtt%qO*X;NGPD`ko)$heFHEojd0WIv^QY z*+C4=;3^Esx}{vhfN$Tf0p9&#>_ZE}!;l!g_PhL@&r=O)P+Pt{Al!?vpUq9)_Eu87 zRpx(gjUO!b-8~RA0WxJGsQxnIjmq6M{m`J4jXh?ebti!vSIeWFfGJUW<%k+wo6Zvy zg#&+_83p1#(SdSKgnqAEJZIIDtJ6N+`G0yce_d2UEJQKj=3NgSL_fs-?5s@%g@*&J z3oN4qtcE@S+h+$TdOzHkK=o~Nbnx7Aq`atMpIJ5D8?J#Fbr!{oe6<_6B{NxTJu$oJ z1A3c~ci+N|Cn{`l2sXgTKya>s~8#sKK%tJ-t84ArfN+mB35 zO{Uq6=md2FlSnE#HDng+x9`npZH35a8Exc@jLLyGfBezf1|jCh0S~yzZ`0P^W@$(M z%ia`yw4^!faaf*JYsqrP3O{R7$!jZ&J`L0$WtR^)=G+qVJa9>ftEi|*o# ze*@NXzO^=E9_M_eDa5zP^HfUKwnub&`^9d2&UIOE(9GOU5_!iEps>lg@+*Vze?CRX z!x_o`42O%A7jl!5sKKD2wIbs#8|wS}2yoj;MiwE|RZR1-s~rHc#<{>zN!|nY|8RFg zBiG&HI%FX4HS-{L(as;1;ShFOUz zOz#qm6@Tt!Ghbv#WZgW%IByPWvE7=37}}X?o?@%r zV#>`*rkRb{e5(&nlPvD?b|iBI)=kojiBbGYS8v)}bUt3g#=Lq=FYI;X4m}4p7X}Ks zdJaiW2q-?uU)ku+GLlZyOs1(uAv7^MBXzM5&qpNkQ-KQ^MB=N}L?dbVPaeT|SNwoZFIY%8*pTkH zSRh3#p6#ZhpNse?RGEH&Yo5lA`yWMzmA>Q!^{<$kN8~2MvDU*HV-}hc)~v(lQn|gC zqGZ)&N(6pmr#=O41Z`XkR31P+&|B`pS~RdgIbnR%eJ)?aJvYV*xbi&EK4@ag_N!Pb zRu4>KV4R}*#ehOv^6U``RZ+=LV!*5S-*-E$;Xgb2_NH(pCirF2QQa_D;mmChS0og6=L7c@Cn%YiW zX6;ABxz$5=>iqh+8D~iTaO>+7)((4D5yXk)^lIQ=J>Wq~!{y>b7r408Kp>*0oQ4j> zS$Qwd$mn#gsr~0!E5i@6HR|776)DAg+IxX>AM=rr7cEw9*x|YA!y6yFXIVmCC*68` zEbKSpruIR8L*(%YF9~QI>&=b00@dOo-m4c$kBu(3bw7J@ukbn9+dj3M?bm1tBnIR- zcFT%we@SCvx}-oU`RavpZ3hQ}ZAM>7;K=Z)?!%MCcCI~E*t*_IYR4-1c71()t*eJo zDBk8A<&a>F3ccHd@8;E$8@aI3gpW zR`;4lw&U^12__R}KhYIGV*sq!$kJ{Zb~~gR41MAlD#1th*Q?;Ol_(eQz5pfNYL2mBy`*ZtXA%5@&xd`J-fB-SXZ8P; z&x&=N_~e9!YTJJn+A$L6pTRBK+(x31cJb8-Jy zTfk)PHuVSd!t-a>sGd<10Pw7xGvr%UIb7S@d)_PiH$E6L%&rrR=Wg;oANFV4C{suS z{|tFu@7K%oRV(yB$pIdP_tKJB#c`L+!P>0aLB%vBaC%MAG=tp)UC3`Zu*vzN*;i9y z_{&^`YY|D0w}tcb?ec6!g)&p_;UE7r+%d0b9ysVU<}fTLVpca$Sf0CnVQ&BWd)H*g zu3RIBYD5$?$20Z;Xlzh0tb7>FFK8|7C#mzduI&kTGSC--#dEBh`n%HZ(1()$N^jUu zJDg@dsC+6|{N*^cVZZv&^HD-s&ZpmX{RPL@TW5XYO;{Fh-AA%_T`FC_OtUI*^SnS*9KWf4Dvc-T2~xWy!;QEFd7gXsT|l^^`D4)f>H53pfdEt?Ug=Ta~Pco%22C&@k@vE_ldwainaZYWL3cEa`{N zng~LQX_MIPYqxh{Q}&hHkGOrlySuxT$e8SD&vOMadx5e6t0a%Lx4J1Ee3K9gH3DBw z{`lwnOKAm$+s{_DtJiF7XLC#%|4O<~C+e2eUT00#%XTy3)41FS5~-IdO`i1HS>lQp zvij61Bp!a+=`60l#8hZ_>;9qo!Lz@gHCK*3N1x^?wQ6jSP+jZQ5%c);FYLQo602JZ zkCmc2PjT(_JmUgVm5K&|t7gTW*svqiS^OjA)>maMr@GvAUscTY|YSL%Jg@{S@W{I#+N%@87@1`We>NL z#N8S#-=Zn+4qn#IPI3{Rzt^?-A!u}RjoA74J-4CqmHMo0L6zZ;8Xk=yahf4H|!2 zOzEYjo}hH;w_312{$`!{1r)~Je@&Gm+nj)<=%fDWLes>d{q>oY*+}u*C2B{{fo>td+w{|ekxze+<_j6cboD7cb)`h{)xw-30 zBVOxOJeBlDE)?f{gGd!Tm^647nS!M@0u2LDQ>FtZP241I_?O|Ejn?3n= z_<=_;0lh z&lMk^Kf9VZDC3X!4%Qa;;{~S73m5;1NF?l`bL3rp^*YJIy0VK@efto_IaKjN0pDV> zv)bulYhQg(FmA&x9&4#RGjS_Pl~ru9ipk5)Pv+;G>Qr(){ry!0##blo<5#z6RNv%Z zP48roWg@8ZxmhaYrVl-6%1Mzoky_b=Cv7wcbxc$t}*0HHzW zR{etw({sj4olaDgfKTEOOT++es5moFX@!8bT=IGi`5nn?X0gRSKfH1@=7r9W1lwE^ zAsX~miOV$}jj+f&t~qR-m*D$KTqo1Mt-`3p&c}`ooa&GCmSM8tz4o;e9n!eT!~ zIf9lbnLkF4Ae5_$gE#YW#|6&rBW&Th2{T%o{Bi3mp;iS}#a7i;ci=8*+HbZq3%7t; zcMC2m4a~jYxW=DVp^@wUlG*K6#_H7Py*Dbm&fkj%1%`{qPP!+|d$Bg34qOi9w*J<) zgxEj^RA2?SRGB_G(9aS&@xyBN9jS;3Z&VR99y^L|Wb)yEb{7v^W&Cp^0t?ZyFB_Ev zPg$5fhom&^K-6T|@~zddH6+j{@nmE^rK*rGs7nA`xbw1gd$06Hc`sHjp;#a{iT`;90xj6yn@_fm;ZCaT(a*{WOat)hsp{+BH7Q26m9IHw300ZJl(;UB z@~%|7WKO?K{_l?{b&|2cEzD-14@44s{?FKPF2_E8e4r*Ov}7T#>Cc|rd80MTCz~!! z=$(IVJ9d>a`Nw~KI`WKz*=sRS>>CUv4U`P;ABzvG59~~K%sLK76Mx|!A~;^@d`^kH z8~3a;v$>lMuOk9Jg0lzXsMqw|bBCq9?;a&O|7{tZZF!f2!)C*l!=A%AKcLT{VBBPkrs=u6i&g5#j)59q&i2u4>8MyAU(~jBZ)Y>> z|2x&((ILO@1ZJ*DG4Db^CYUV05wF>)(c7>8qnZz9VzBtI$SN34|L~K>5I{?AON|M+ zG6Vjxqc{c7f|?(mBAVbaaSFrH(S=MfIIzUUorBhyp8@~mOtXFFwx&e~N{(Om zfMB8-_bY1G%Wa$Bo(@@)?T#D!uDPERQEl;!?oa-o#_2BN(k1))7ChvT*9d;cp`XFx zaC=CKztH!pK?4_^)8e$QBdr|tgMU)6`$m?&s7kdr%i`t1rAqmr39X?|Xw#)cfKM+yE&qiIa^ZCct4LrYj+ zyZ*cV%26LwMhpt=g9(@WzHHUOg3PtB?|%yBpEP|%l_$W@>9voar9EHz#n2K}&uu!~ zSIf}yI%E!IfP5nn^Y^OegF5QLg8c~!i>mmjLs_L)>|Dth4*bxE;ai4@?dImy9kHBY ziDZ(mj*dn{m&Rfm;-KYs&~kQP-?N#=_H^XrT83VmyJANJRSYdYQ6Cvv4h}T6wXJ)E zXi<)rXi$Z{I|=r?^RpIYLz}lppYPAK&pL-c(9()YNbrC23FE^vbWHE!GBN(R2&Rto zuE*2=u&w_b`Bb$#y9GL@rq5O}90gtc*^Y^WvdXLS6utSDSkI{Xoh)`zYqWak!J>ze zk^0us=zBSDlmjvw_uUU#cJ<4};@*i;wXHcDl`UQ)?Ls4@LH#Z;%HjG`|7-RY!&Lgk zS)W>;Z#kmZ*A*|Jp`WzAo%;E6MQ@^_#j>`7kt>tB gi1Yu$AOE3l7D1bvv-j14Z$K%@tI6faJb3ZH07pZiU;qFB literal 0 HcmV?d00001 diff --git a/docs/visualisation.rst b/docs/visualisation.rst index e7a03336..45e5763d 100644 --- a/docs/visualisation.rst +++ b/docs/visualisation.rst @@ -4,6 +4,14 @@ It is quite often that the user want to visualise the results to gain confidence on the computed free energy. **alchemlyb** provides various visualisation tools to help user to judge the estimate. +.. currentmodule:: alchemlyb.visualisation + +.. autosummary:: + :toctree: visualisation + + plot_mbar_overlap_matrix + plot_ti_dhdl + .. _plot_overlap_matrix: Overlap Matrix of the MBAR @@ -34,6 +42,38 @@ Will give a plot looks like this .. image:: images/O_MBAR.png +.. _plot_TI_dhdl: + +dhdl Plot of the TI +------------------- +In order for the :class:`~alchemlyb.estimators.TI` estimator to work reliably, +the change in the dhdl between lambda state 0 and lambda state 1 should be +adequately sampled. The function :func:`~alchemlyb.visualisation.plot_ti_dhdl` +can be used to access the change of the dhdl across the lambda states. + +More than one :class:`~alchemlyb.estimators.TI` estimators can be plotted +together as well. :: + + >>> import pandas as pd + >>> from alchemtest.gmx import load_benzene + >>> from alchemlyb.parsing.gmx import extract_dHdl + >>> from alchemlyb.estimators import TI + + >>> bz = load_benzene().data + >>> dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']]) + >>> ti_coul = TI().fit(dHdl_coul) + >>> dHdl_vdw = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']]) + >>> ti_vdw = TI().fit(dHdl_vdw) + + >>> from alchemlyb.visualisation import plot_ti_dhdl + >>> ax = plot_ti_dhdl([ti_coul, ti_vdw], labels=['Coul', 'VDW'], colors=['r', 'g']) + >>> ax.figure.savefig('dhdl_TI.pdf') + + +Will give a plot looks like this + +.. image:: images/dhdl_TI.png + .. [Klimovich2015] Klimovich, P.V., Shirts, M.R. & Mobley, D.L. Guidelines for the analysis of free energy calculations. J Comput Aided Mol Des 29, 397–411 (2015). https://doi.org/10.1007/s10822-015-9840-9 diff --git a/docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst b/docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst index 5323d1a9..a4de37da 100644 --- a/docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst +++ b/docs/visualisation/alchemlyb.visualisation.plot_mbar_overlap_matrix.rst @@ -1,4 +1,4 @@ -.. _plot_overlap: +.. _visualisation_plot_mbar_overlap_matrix: Plot Overlap Matrix from MBAR ============================= diff --git a/docs/visualisation/alchemlyb.visualisation.plot_ti_dhdl.rst b/docs/visualisation/alchemlyb.visualisation.plot_ti_dhdl.rst new file mode 100644 index 00000000..d4247515 --- /dev/null +++ b/docs/visualisation/alchemlyb.visualisation.plot_ti_dhdl.rst @@ -0,0 +1,22 @@ +.. _visualisation_plot_mbar_plot_ti_dhdl: + +Plot dhdl from TI +================= + +The function :func:`~alchemlyb.visualisation.plot_ti_dhdl` allows +the user to plot the dhdl from :class:`~alchemlyb.estimators.TI` estimator. +Several :class:`~alchemlyb.estimators.TI` estimators could be passed to the +function to give a concerted picture of the whole alchemical transformation. +When custom labels are desirable, the user could pass a list of strings to the +*labels* for labelling each alchemical transformation differently. The color of +each alchemical transformation could also be set by passing a list of color +string to the *colors*. The unit in the y axis could be labelled to other units +by setting *units*, which by default is kcal/mol. The user can pass +:class:`matplotlib.axes.Axes` into the function to have the dhdl drawn on a +specific axes. + +Please check :ref:`How to plot TI dhdl ` for usage. + +API Reference +------------- +.. autofunction:: alchemlyb.visualisation.plot_ti_dhdl \ No newline at end of file diff --git a/src/alchemlyb/estimators/mbar_.py b/src/alchemlyb/estimators/mbar_.py index dd8127ae..975d177b 100644 --- a/src/alchemlyb/estimators/mbar_.py +++ b/src/alchemlyb/estimators/mbar_.py @@ -45,9 +45,6 @@ class MBAR(BaseEstimator): states_ : list Lambda states for which free energy differences were obtained. - overlap_matrix: numpy.matrix - The overlap matrix. - """ def __init__(self, maximum_iterations=10000, relative_tolerance=1.0e-7, diff --git a/src/alchemlyb/estimators/ti_.py b/src/alchemlyb/estimators/ti_.py index 7f25dbe3..383341c9 100644 --- a/src/alchemlyb/estimators/ti_.py +++ b/src/alchemlyb/estimators/ti_.py @@ -26,6 +26,9 @@ class TI(BaseEstimator): states_ : list Lambda states for which free energy differences were obtained. + dhdl : DataFrame + The estimated dhdl of each state. + """ def __init__(self, verbose=False): @@ -92,6 +95,7 @@ def fit(self, dHdl): self.delta_f_ = pd.DataFrame(adelta - adelta.T, columns=means.index.values, index=means.index.values) + self.dhdl = means # yield standard deviation d_delta_f_ between each state self.d_delta_f_ = pd.DataFrame(np.sqrt(ad_delta + ad_delta.T), diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index d482caad..e9b5434c 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -2,9 +2,10 @@ import pandas as pd from alchemtest.gmx import load_benzene -from alchemlyb.parsing.gmx import extract_u_nk -from alchemlyb.estimators import MBAR +from alchemlyb.parsing.gmx import extract_u_nk, extract_dHdl +from alchemlyb.estimators import MBAR, TI from alchemlyb.visualisation.mbar_matrix import plot_mbar_overlap_matrix +from alchemlyb.visualisation.ti_dhdl import plot_ti_dhdl def test_plot_mbar_omatrix(): '''Just test if the plot runs''' @@ -25,3 +26,20 @@ def test_plot_mbar_omatrix(): assert isinstance(plot_mbar_overlap_matrix(overlap_maxtrix), matplotlib.axes.Axes) +def test_plot_ti_dhdl(): + '''Just test if the plot runs''' + bz = load_benzene().data + dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']]) + ti_coul = TI() + ti_coul.fit(dHdl_coul) + assert isinstance(plot_ti_dhdl(ti_coul), + matplotlib.axes.Axes) + assert isinstance(plot_ti_dhdl(ti_coul, labels=['Coul']), + matplotlib.axes.Axes) + assert isinstance(plot_ti_dhdl(ti_coul, labels=['Coul'], colors=['r']), + matplotlib.axes.Axes) + dHdl_vdw = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']]) + ti_vdw = TI().fit(dHdl_vdw) + assert isinstance(plot_ti_dhdl([ti_coul, ti_vdw]), + matplotlib.axes.Axes) + diff --git a/src/alchemlyb/visualisation/__init__.py b/src/alchemlyb/visualisation/__init__.py index 4fc73051..5cf2c3f9 100644 --- a/src/alchemlyb/visualisation/__init__.py +++ b/src/alchemlyb/visualisation/__init__.py @@ -1 +1,2 @@ -from .mbar_matrix import plot_mbar_overlap_matrix \ No newline at end of file +from .mbar_matrix import plot_mbar_overlap_matrix +from .ti_dhdl import plot_ti_dhdl \ No newline at end of file diff --git a/src/alchemlyb/visualisation/ti_dhdl.py b/src/alchemlyb/visualisation/ti_dhdl.py new file mode 100644 index 00000000..916455db --- /dev/null +++ b/src/alchemlyb/visualisation/ti_dhdl.py @@ -0,0 +1,169 @@ +"""Functions for Plotting the dhdl for the TI estimator. + +To assess the quality of the TI estimator, the dhdl from lambda state 0 +to lambda state 1 can plotted to assess if the change in dhdl is sampled +thoroughly. + +The code for producing the dhdl plot is modified based on +: `Alchemical Analysis `_. + +""" +from __future__ import division + +import matplotlib.pyplot as plt +from matplotlib.font_manager import FontProperties as FP +import numpy as np + +def plot_ti_dhdl(dhdl_data, labels=None, colors=None, units='kcal/mol', ax=None): + '''Plot the dhdl of TI. + + Parameters + ---------- + dhdl_data : :class:`~alchemlyb.estimators.TI` or list + One or more :class:`~alchemlyb.estimators.TI` estimator, where the + dhdl value will be taken from. + labels : List + list of labels for labelling all the alchemical transformations. + colors : List + list of colors for plotting all the alchemical transformations. + Default: ['r', 'g', '#7F38EC', '#9F000F', 'b', 'y'] + units : str + The unit of the estimate. Default: 'kcal/mol' + ax : matplotlib.axes.Axes + Matplotlib axes object where the plot will be drawn on. If ax=None, + a new axes will be generated. + + Returns + ------- + matplotlib.axes.Axes + An axes with the overlap matrix drawn. + + Notes + ----- + The code is taken and modified from + : `Alchemical Analysis `_ + + ''' + # Make it into a list + try: + len(dhdl_data) + except TypeError: + dhdl_data = [dhdl_data, ] + + if ax is None: + fig, ax = plt.subplots(figsize=(8, 6)) + + ax.spines['bottom'].set_position('zero') + ax.spines['top'].set_color('none') + ax.spines['right'].set_color('none') + ax.xaxis.set_ticks_position('bottom') + ax.yaxis.set_ticks_position('left') + + for k, spine in ax.spines.items(): + spine.set_zorder(12.2) + + # Make level names + if labels is None: + lv_names2 = [] + for dhdl in dhdl_data: + # Assume that the dhdl has only one columns + lv_names2.append(dhdl.dhdl.columns.values[0].capitalize()) + else: + assert len(labels) == len(dhdl_data), \ + 'Length of labels ({}) should be the same as the number of data ({})'.format(len(labels), len(dhdl_data)) + lv_names2 = labels + + if colors is None: + colors = ['r', 'g', '#7F38EC', '#9F000F', 'b', 'y'] + else: + assert len(colors) >= len(dhdl_data), \ + 'Number of colors ({}) should be larger than the number of data ({})'.format(len(labels), len(dhdl_data)) + # Get the real data out + xs, ndx, dx = [0], 0, 0.001 + min_y, max_y = 0, 0 + K = -1 + for dhdl in dhdl_data: + x = dhdl.dhdl.index.values + K += len(x) + y = dhdl.dhdl.values.ravel() + + min_y = min(y.min(), min_y) + max_y = max(y.max(), max_y) + + for i in range(len(x) - 1): + # pl.plot(x,y) + if i % 2 == 0: + ax.fill_between(x[i:i + 2] + ndx, 0, y[i:i + 2], + color=colors[ndx], alpha=1.0) + else: + ax.fill_between(x[i:i + 2] + ndx, 0, y[i:i + 2], + color=colors[ndx], alpha=0.5) + + xlegend = [-100 * wnum for wnum in range(len(lv_names2))] + ax.plot(xlegend, [0 * wnum for wnum in xlegend], ls='-', + color=colors[ndx], + label=lv_names2[ndx]) ## for the paper + xs += (x + ndx).tolist()[1:] + ndx += 1 + # + # + # min_dl = dlam[dlam != 0].min() + # S = int(0.4 / min_dl) + # + # Make sure the tick labels are not overcrowded. + xs = np.array(xs) + dl_mat = np.array([xs - i for i in xs]) + ri = range(len(xs)) + + def getInd(r=ri, z=[0]): + primo = r[0] + min_dl = ndx * 0.02 * 2 ** (primo > 10) + if dl_mat[primo].max() < min_dl: + return z + for i in r: + for j in range(len(xs)): + if dl_mat[i, j] > min_dl: + z.append(j) + return getInd(ri[j:], z) + + xt = [i if (i in getInd()) else '' for i in range(K)] + ax.set_xticks(xs[1:],) + ax.set_xticklabels(xt[1:], fontsize=10) + ax.yaxis.label.set_size(10) + + # Remove the abscissa ticks and set up the axes limits. + for tick in ax.get_xticklines(): + tick.set_visible(False) + ax.set_xlim(0, ndx) + min_y *= 1.01 + max_y *= 1.01 + + # Modified so that the x label won't conflict with the lambda label + min_y -= (max_y-min_y)*0.1 + + ax.set_ylim(min_y, max_y) + + for i, j in zip(xs[1:], xt[1:]): + ax.annotate( + ('%.2f' % (i - 1.0 if i > 1.0 else i) if not j == '' else ''), + xy=(i, 0), xytext=(i, 0.01), size=10, rotation=90, + textcoords=('data', 'axes fraction'), va='bottom', ha='center', + color='#151B54') + if ndx > 1: + lenticks = len(ax.get_ymajorticklabels()) - 1 + if min_y < 0: lenticks -= 1 + if lenticks < 5: + from matplotlib.ticker import AutoMinorLocator as AML + ax.yaxis.set_minor_locator(AML()) + ax.grid(which='both', color='w', lw=0.25, axis='y', zorder=12) + ax.set_ylabel( + r'$\mathrm{\langle{\frac{ \partial U } { \partial \lambda }}\rangle_{\lambda}\/%s}$' % units, + fontsize=20, color='#151B54') + ax.annotate('$\mathit{\lambda}$', xy=(0, 0), xytext=(0.5, -0.05), size=18, + textcoords='axes fraction', va='top', ha='center', + color='#151B54') + lege = ax.legend(prop=FP(size=14), frameon=False, loc=1) + for l in lege.legendHandles: + l.set_linewidth(10) + return ax + From 3efbaefbfbd3fe7c2d4935e640ac7f155c62e4ce Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 21:05:03 +0100 Subject: [PATCH 18/25] fix test --- src/alchemlyb/visualisation/ti_dhdl.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/alchemlyb/visualisation/ti_dhdl.py b/src/alchemlyb/visualisation/ti_dhdl.py index 916455db..aa324e9a 100644 --- a/src/alchemlyb/visualisation/ti_dhdl.py +++ b/src/alchemlyb/visualisation/ti_dhdl.py @@ -81,10 +81,8 @@ def plot_ti_dhdl(dhdl_data, labels=None, colors=None, units='kcal/mol', ax=None) # Get the real data out xs, ndx, dx = [0], 0, 0.001 min_y, max_y = 0, 0 - K = -1 for dhdl in dhdl_data: x = dhdl.dhdl.index.values - K += len(x) y = dhdl.dhdl.values.ravel() min_y = min(y.min(), min_y) @@ -126,9 +124,14 @@ def getInd(r=ri, z=[0]): z.append(j) return getInd(ri[j:], z) - xt = [i if (i in getInd()) else '' for i in range(K)] - ax.set_xticks(xs[1:],) - ax.set_xticklabels(xt[1:], fontsize=10) + xt = [] + for i in range(len(xs)): + if i in getInd(): + xt.append(i) + else: + xt.append('') + + plt.xticks(xs[1:], xt[1:], fontsize=10) ax.yaxis.label.set_size(10) # Remove the abscissa ticks and set up the axes limits. From d68ba2f390b58dba343dc704648f65c91c43868b Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 21:12:06 +0100 Subject: [PATCH 19/25] bump test --- src/alchemlyb/tests/test_visualisation.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index e9b5434c..90389c78 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -1,4 +1,5 @@ import matplotlib +import matplotlib.pyplot as plt import pandas as pd from alchemtest.gmx import load_benzene @@ -34,6 +35,9 @@ def test_plot_ti_dhdl(): ti_coul.fit(dHdl_coul) assert isinstance(plot_ti_dhdl(ti_coul), matplotlib.axes.Axes) + fig, ax = plt.subplots(figsize=(8, 6)) + assert isinstance(plot_ti_dhdl(ti_coul, ax=ax), + matplotlib.axes.Axes) assert isinstance(plot_ti_dhdl(ti_coul, labels=['Coul']), matplotlib.axes.Axes) assert isinstance(plot_ti_dhdl(ti_coul, labels=['Coul'], colors=['r']), @@ -42,4 +46,9 @@ def test_plot_ti_dhdl(): ti_vdw = TI().fit(dHdl_vdw) assert isinstance(plot_ti_dhdl([ti_coul, ti_vdw]), matplotlib.axes.Axes) + dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb'][:3]]) + ti_coul = TI() + ti_coul.fit(dHdl_coul) + assert isinstance(plot_ti_dhdl(ti_coul), + matplotlib.axes.Axes) From 1456ac2ab84a7c4e59a8e7fce60d98a7798ba832 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 21:23:46 +0100 Subject: [PATCH 20/25] bump test --- src/alchemlyb/tests/test_visualisation.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/alchemlyb/tests/test_visualisation.py b/src/alchemlyb/tests/test_visualisation.py index 90389c78..06808ea6 100644 --- a/src/alchemlyb/tests/test_visualisation.py +++ b/src/alchemlyb/tests/test_visualisation.py @@ -1,6 +1,7 @@ import matplotlib import matplotlib.pyplot as plt import pandas as pd +import numpy as np from alchemtest.gmx import load_benzene from alchemlyb.parsing.gmx import extract_u_nk, extract_dHdl @@ -46,9 +47,10 @@ def test_plot_ti_dhdl(): ti_vdw = TI().fit(dHdl_vdw) assert isinstance(plot_ti_dhdl([ti_coul, ti_vdw]), matplotlib.axes.Axes) - dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb'][:3]]) - ti_coul = TI() - ti_coul.fit(dHdl_coul) + ti_coul.dhdl = pd.DataFrame.from_dict( + {'fep': range(100)}, + orient='index', + columns=np.arange(100)/100).T assert isinstance(plot_ti_dhdl(ti_coul), matplotlib.axes.Axes) From 8de28f5811999617c7a97c9819b41cfe341eaf24 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 22:40:57 +0100 Subject: [PATCH 21/25] # pragma: no cover --- src/alchemlyb/visualisation/ti_dhdl.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/alchemlyb/visualisation/ti_dhdl.py b/src/alchemlyb/visualisation/ti_dhdl.py index aa324e9a..8c582214 100644 --- a/src/alchemlyb/visualisation/ti_dhdl.py +++ b/src/alchemlyb/visualisation/ti_dhdl.py @@ -118,7 +118,7 @@ def getInd(r=ri, z=[0]): min_dl = ndx * 0.02 * 2 ** (primo > 10) if dl_mat[primo].max() < min_dl: return z - for i in r: + for i in r: # pragma: no cover for j in range(len(xs)): if dl_mat[i, j] > min_dl: z.append(j) @@ -155,7 +155,7 @@ def getInd(r=ri, z=[0]): if ndx > 1: lenticks = len(ax.get_ymajorticklabels()) - 1 if min_y < 0: lenticks -= 1 - if lenticks < 5: + if lenticks < 5: # pragma: no cover from matplotlib.ticker import AutoMinorLocator as AML ax.yaxis.set_minor_locator(AML()) ax.grid(which='both', color='w', lw=0.25, axis='y', zorder=12) From 97962a6aa02886cf7621e6e839f187f189a1ae4f Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 22:45:33 +0100 Subject: [PATCH 22/25] typo --- docs/visualisation.rst | 2 +- src/alchemlyb/visualisation/ti_dhdl.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/visualisation.rst b/docs/visualisation.rst index 45e5763d..8d67d3b3 100644 --- a/docs/visualisation.rst +++ b/docs/visualisation.rst @@ -49,7 +49,7 @@ dhdl Plot of the TI In order for the :class:`~alchemlyb.estimators.TI` estimator to work reliably, the change in the dhdl between lambda state 0 and lambda state 1 should be adequately sampled. The function :func:`~alchemlyb.visualisation.plot_ti_dhdl` -can be used to access the change of the dhdl across the lambda states. +can be used to assess the change of the dhdl across the lambda states. More than one :class:`~alchemlyb.estimators.TI` estimators can be plotted together as well. :: diff --git a/src/alchemlyb/visualisation/ti_dhdl.py b/src/alchemlyb/visualisation/ti_dhdl.py index 8c582214..61d676a6 100644 --- a/src/alchemlyb/visualisation/ti_dhdl.py +++ b/src/alchemlyb/visualisation/ti_dhdl.py @@ -36,7 +36,7 @@ def plot_ti_dhdl(dhdl_data, labels=None, colors=None, units='kcal/mol', ax=None) Returns ------- matplotlib.axes.Axes - An axes with the overlap matrix drawn. + An axes with the TI dhdl drawn. Notes ----- From 3f9a651efea6bc6eb538affa633fc2fc14bc59c4 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Thu, 1 Oct 2020 22:54:52 +0100 Subject: [PATCH 23/25] remove old comments --- src/alchemlyb/visualisation/ti_dhdl.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/alchemlyb/visualisation/ti_dhdl.py b/src/alchemlyb/visualisation/ti_dhdl.py index 61d676a6..bf5494a3 100644 --- a/src/alchemlyb/visualisation/ti_dhdl.py +++ b/src/alchemlyb/visualisation/ti_dhdl.py @@ -89,7 +89,6 @@ def plot_ti_dhdl(dhdl_data, labels=None, colors=None, units='kcal/mol', ax=None) max_y = max(y.max(), max_y) for i in range(len(x) - 1): - # pl.plot(x,y) if i % 2 == 0: ax.fill_between(x[i:i + 2] + ndx, 0, y[i:i + 2], color=colors[ndx], alpha=1.0) @@ -100,14 +99,10 @@ def plot_ti_dhdl(dhdl_data, labels=None, colors=None, units='kcal/mol', ax=None) xlegend = [-100 * wnum for wnum in range(len(lv_names2))] ax.plot(xlegend, [0 * wnum for wnum in xlegend], ls='-', color=colors[ndx], - label=lv_names2[ndx]) ## for the paper + label=lv_names2[ndx]) xs += (x + ndx).tolist()[1:] ndx += 1 - # - # - # min_dl = dlam[dlam != 0].min() - # S = int(0.4 / min_dl) - # + # Make sure the tick labels are not overcrowded. xs = np.array(xs) dl_mat = np.array([xs - i for i in xs]) From 3a2e2356174e9d03a6fed381267dc15bcf37dcab Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Fri, 2 Oct 2020 09:37:00 +0100 Subject: [PATCH 24/25] Update ti_dhdl.py --- src/alchemlyb/visualisation/ti_dhdl.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/alchemlyb/visualisation/ti_dhdl.py b/src/alchemlyb/visualisation/ti_dhdl.py index bf5494a3..7b53217b 100644 --- a/src/alchemlyb/visualisation/ti_dhdl.py +++ b/src/alchemlyb/visualisation/ti_dhdl.py @@ -69,15 +69,23 @@ def plot_ti_dhdl(dhdl_data, labels=None, colors=None, units='kcal/mol', ax=None) # Assume that the dhdl has only one columns lv_names2.append(dhdl.dhdl.columns.values[0].capitalize()) else: - assert len(labels) == len(dhdl_data), \ - 'Length of labels ({}) should be the same as the number of data ({})'.format(len(labels), len(dhdl_data)) - lv_names2 = labels + if len(labels) == len(dhdl_data): + lv_names2 = labels + else: # pragma: no cover + raise ValueError( + 'Length of labels ({}) should be the same as the number of data ({})'.format( + len(labels), len(dhdl_data))) if colors is None: colors = ['r', 'g', '#7F38EC', '#9F000F', 'b', 'y'] else: - assert len(colors) >= len(dhdl_data), \ - 'Number of colors ({}) should be larger than the number of data ({})'.format(len(labels), len(dhdl_data)) + if len(colors) >= len(dhdl_data): + pass + else: # pragma: no cover + raise ValueError( + 'Number of colors ({}) should be larger than the number of data ({})'.format( + len(labels), len(dhdl_data))) + # Get the real data out xs, ndx, dx = [0], 0, 0.001 min_y, max_y = 0, 0 From f7853fcc0e66e3565fd0b6c1d65cd70132471ae7 Mon Sep 17 00:00:00 2001 From: Zhiyi Wu Date: Fri, 2 Oct 2020 09:38:20 +0100 Subject: [PATCH 25/25] Update CHANGES --- CHANGES | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGES b/CHANGES index b95aa851..c64a912b 100644 --- a/CHANGES +++ b/CHANGES @@ -21,6 +21,7 @@ The rules for this file: Enhancements - Allow the overlap matrix of the MBAR estimator to be plotted. (PR #107) + - Allow the dhdl of the TI estimator to be plotted. (PR #110) Deprecations