Allow dF-state to be plotted (#112)

* partially address #73 
* add df_state plotting code in in new visualisation.dF_state, based on alchemical-analysis.py
* add test
* Update CHANGES

CHANGES

@@ -20,8 +20,9 @@ The rules for this file:
 * 0.?.?
 
 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)
+  - Allow the overlap matrix of the MBAR estimator to be plotted. (issue #73, PR #107)
+  - Allow the dhdl of the TI estimator to be plotted. (issue #73, PR #110)
+  - Allow the dF states to be plotted. (issue #73, PR #112)
 
 Deprecations
 

docs/visualisation.rst

@@ -11,6 +11,7 @@ visualisation tools to help user to judge the estimate.
 
     plot_mbar_overlap_matrix
     plot_ti_dhdl
+    plot_dF_state
 
 .. _plot_overlap_matrix:
 
@@ -40,7 +41,12 @@ the degree of overlap. It is recommended that there should be at least
 
 Will give a plot looks like this
 
-.. image:: images/O_MBAR.png
+.. figure:: images/O_MBAR.png
+
+   Overlap between the distributions of potential energy differences is
+   essential for accurate free energy calculations and can be quantified by
+   computing the overlap matrix 𝐎. Its elements 𝑂𝑖𝑗 are the probabilities of
+   observing a sample from state i (𝑖 th row) in state j (𝑗 th column).
 
 .. _plot_TI_dhdl:
 
@@ -72,7 +78,59 @@ together as well. ::
 
 Will give a plot looks like this
 
-.. image:: images/dhdl_TI.png
+.. figure:: images/dhdl_TI.png
+
+   A plot of ⟨∂𝑈/∂𝜆⟩ versus 𝜆 for thermodynamic integration, with filled areas
+   indicating free energy estimates from the trapezoid rule. Different 𝛥𝐺
+   components are shown in distinct colors: in red is the electrostatic 𝛥𝐺
+   component (𝜆 indices 0–4), while in green is the van der Waals 𝛥𝐺 component
+   (𝜆 indices 5–19). Color intensity alternates with increasing 𝜆 index.
+
+.. _plot_dF_states:
+
+dF States Plots between Different estimators
+--------------------------------------------
+Another way of assessing the quality of free energy estimate would be comparing
+the free energy difference between adjacent lambda states (dF) using different
+estimators [Klimovich2015]_. The function :func:`~alchemlyb.visualisation.plot_dF_state` can
+be used, for example, to compare the dF of both Coulombic and VDW
+transformations using :class:`~alchemlyb.estimators.TI`,
+:class:`~alchemlyb.estimators.BAR` and :class:`~alchemlyb.estimators.MBAR`
+estimators. ::
+
+    >>> from alchemtest.gmx import load_benzene
+    >>> from alchemlyb.parsing.gmx import extract_u_nk, extract_dHdl
+    >>> from alchemlyb.estimators import MBAR, TI, BAR
+    >>> import matplotlib.pyplot as plt
+    >>> import pandas as pd
+    >>> from alchemlyb.visualisation.dF_state import plot_dF_state
+    >>> bz = load_benzene().data
+    >>> u_nk_coul = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']])
+    >>> dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']])
+    >>> u_nk_vdw = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['VDW']])
+    >>> dHdl_vdw = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']])
+    >>> ti_coul = TI().fit(dHdl_coul)
+    >>> ti_vdw = TI().fit(dHdl_vdw)
+    >>> bar_coul = BAR().fit(u_nk_coul)
+    >>> bar_vdw = BAR().fit(u_nk_vdw)
+    >>> mbar_coul = MBAR().fit(u_nk_coul)
+    >>> mbar_vdw = MBAR().fit(u_nk_vdw)
+
+    >>> estimators = [(ti_coul, ti_vdw),
+                     (bar_coul, bar_vdw),
+                     (mbar_coul, mbar_vdw),]
+
+    >>> fig = plot_dF_state(estimators, orientation='portrait')
+    >>> fig.savefig('dF_state.pdf', bbox_inches='tight')
+
+
+Will give a plot looks like this
+
+.. figure:: images/dF_states.png
+
+   A bar plot of the free energy differences evaluated between pairs of adjacent
+   states via several methods, with corresponding error estimates for each method.
+
 
 .. [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

docs/visualisation/alchemlyb.visualisation.plot_dF_state.rst

@@ -0,0 +1,35 @@
+.. _visualisation_plot_dF_state:
+
+Plot dF states from multiple estimators
+=======================================
+
+The function :func:`~alchemlyb.visualisation.plot_dF_state` allows the user to
+plot and compare the free energy difference between states ("dF") from various
+kinds of :class:`~alchemlyb.estimators`.
+
+To compare the dF states of a single alchemical transformation among various
+:class:`~alchemlyb.estimators`, the user can pass a list of `estimators`. (e.g.
+`estimators` = [:class:`~alchemlyb.estimators.TI`,
+:class:`~alchemlyb.estimators.BAR`, :class:`~alchemlyb.estimators.MBAR`])
+
+To compare the dF states of a multiple alchemical transformations, results from
+the same :class:`~alchemlyb.estimators` can be concatenated into a list, which
+is then bundled to to another list of different :class:`~alchemlyb.estimators`.
+(e.g. `estimators` = [(:class:`~alchemlyb.estimators.TI`,
+:class:`~alchemlyb.estimators.TI`), (:class:`~alchemlyb.estimators.BAR`,
+:class:`~alchemlyb.estimators.BAR`), (:class:`~alchemlyb.estimators.MBAR`,
+:class:`~alchemlyb.estimators.MBAR`)])
+
+The figure could be plotted in *portrait* or *landscape* mode by setting the
+`orientation`. `nb` is used to control the number of dF states in one row.
+The user could pass a list of strings to `labels` to name the
+:class:`~alchemlyb.estimators` or a list of strings to `colors` to color
+the estimators differently. The unit in the y axis could be labelled to other
+units by setting `units`, which by default is kcal/mol.
+
+Please check :ref:`How to plot dF states <plot_dF_states>` for a complete
+example.
+
+API Reference
+-------------
+.. autofunction:: alchemlyb.visualisation.plot_dF_state

src/alchemlyb/tests/test_visualisation.py

@@ -2,12 +2,14 @@
 import matplotlib.pyplot as plt
 import pandas as pd
 import numpy as np
+import pytest
 
 from alchemtest.gmx import load_benzene
 from alchemlyb.parsing.gmx import extract_u_nk, extract_dHdl
-from alchemlyb.estimators import MBAR, TI
+from alchemlyb.estimators import MBAR, TI, BAR
 from alchemlyb.visualisation.mbar_matrix import plot_mbar_overlap_matrix
 from alchemlyb.visualisation.ti_dhdl import plot_ti_dhdl
+from alchemlyb.visualisation.dF_state import plot_dF_state
 
 def test_plot_mbar_omatrix():
     '''Just test if the plot runs'''
@@ -54,3 +56,43 @@ def test_plot_ti_dhdl():
     assert isinstance(plot_ti_dhdl(ti_coul),
                matplotlib.axes.Axes)
 
+def test_plot_dF_state():
+    '''Just test if the plot runs'''
+    bz = load_benzene().data
+    u_nk_coul = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['Coulomb']])
+    dHdl_coul = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['Coulomb']])
+    u_nk_vdw = pd.concat([extract_u_nk(xvg, T=300) for xvg in bz['VDW']])
+    dHdl_vdw = pd.concat([extract_dHdl(xvg, T=300) for xvg in bz['VDW']])
+
+    ti_coul = TI().fit(dHdl_coul)
+    ti_vdw = TI().fit(dHdl_vdw)
+    bar_coul = BAR().fit(u_nk_coul)
+    bar_vdw = BAR().fit(u_nk_vdw)
+    mbar_coul = MBAR().fit(u_nk_coul)
+    mbar_vdw = MBAR().fit(u_nk_vdw)
+
+    dhdl_data = [(ti_coul, ti_vdw),
+                 (bar_coul, bar_vdw),
+                 (mbar_coul, mbar_vdw), ]
+    fig = plot_dF_state(dhdl_data, orientation='portrait')
+    assert isinstance(fig, matplotlib.figure.Figure)
+    fig = plot_dF_state(dhdl_data, orientation='landscape')
+    assert isinstance(fig, matplotlib.figure.Figure)
+    fig = plot_dF_state(dhdl_data, labels=['MBAR', 'TI', 'BAR'])
+    assert isinstance(fig, matplotlib.figure.Figure)
+    with pytest.raises(ValueError):
+        fig = plot_dF_state(dhdl_data, labels=['MBAR', 'TI',])
+    fig = plot_dF_state(dhdl_data, colors=['#C45AEC', '#33CC33', '#F87431'])
+    assert isinstance(fig, matplotlib.figure.Figure)
+    with pytest.raises(ValueError):
+        fig = plot_dF_state(dhdl_data, colors=['#C45AEC', '#33CC33'])
+    with pytest.raises(NameError):
+        fig = plot_dF_state(dhdl_data, orientation='xxx')
+    fig = plot_dF_state(ti_coul, orientation='landscape')
+    assert isinstance(fig, matplotlib.figure.Figure)
+    fig = plot_dF_state(ti_coul, orientation='portrait')
+    assert isinstance(fig, matplotlib.figure.Figure)
+    fig = plot_dF_state([ti_coul, bar_coul])
+    assert isinstance(fig, matplotlib.figure.Figure)
+    fig = plot_dF_state([(ti_coul, ti_vdw)])
+    assert isinstance(fig, matplotlib.figure.Figure)

src/alchemlyb/visualisation/__init__.py

@@ -1,2 +1,3 @@
 from .mbar_matrix import plot_mbar_overlap_matrix
-from .ti_dhdl import plot_ti_dhdl
+from .ti_dhdl import plot_ti_dhdl
+from .dF_state import plot_dF_state