Implement the workflow pipeline

CHANGES

@@ -13,14 +13,19 @@ The rules for this file:
   * release numbers follow "Semantic Versioning" https://semver.org
 
 ------------------------------------------------------------------------------
-??/??/2022 orbeckst
+??/??/2022 orbeckst, xiki-tempula
 
   * 1.0.0
 
 Changes
   - Drop support for py3.7 (Issue #179, PR #214)
+  - forward_backward_convergence will use AutoMBAR as backend when `MBAR` is
+    selected as the estimator (PR #114).
+  - AutoMBAR accepts the `method` argument (PR #114).
 
 Enhancements
+  - Add a base class for workflows (PR #188).
+  - Add the ABFE workflow (PR #114).
 
 
 Fixes
@@ -39,7 +44,6 @@ Changes
   - remove broken .zip support from util.anyopen() (PR #197)
 
 Enhancements
-  - Add a base class for workflows (PR #188).
   - Add filter function to gmx.extract to make it more robust (PR #183): can filter 
     incomplete/corrupted lines (#126, #171) with filter=True.
   - Add support to util.anyopen() for taking filelike objects (PR #197)

docs/api_principles.rst

@@ -65,7 +65,8 @@ The library is structured as follows, following a similar style to
     │   ├── ti_dhdl.py
     │   └── ...
     └── workflows           ### WORK IN PROGRESS
-        └── ...
+    │   ├── base.py
+    │   ├── abfe.py
 
 
 * The :mod:`~alchemlyb.parsing` submodule contains parsers for individual MD engines, since the output files needed to perform alchemical free energy calculations vary widely and are not standardized.  Each module at the very least provides an `extract_u_nk` function for extracting reduced potentials (needed for MBAR), as well as an `extract_dHdl` function for extracting derivatives required for thermodynamic integration.  Other helper functions may be exposed for additional processing, such as generating an XVG file from an EDR file in the case of GROMACS.  All `extract\_*` functions take similar arguments (a file path, parameters such as temperature), and produce standard outputs (:class:`pandas.DataFrame` for reduced potentials, :class:`pandas.Series` for derivatives).
@@ -74,8 +75,8 @@ The library is structured as follows, following a similar style to
 * The :mod:`~alchemlyb.convergence` submodule features convenience functions/classes for doing convergence analysis using a given dataset and a chosen estimator.
 * The :mod:`~alchemlyb.postprocessors` submodule contains functions to calculate new quantities or express data in different units.
 * The :mod:`~alchemlyb.visualisation` submodule contains convenience plotting functions as known from, for example, `alchemical-analysis.py`_.
-* The :mod:`~alchemlyb.workflows` submodule will contain complete analysis workflows that will serve as larger building blocks for complex analysis pipelines or a command line interface.
-
+* The :mod:`~alchemlyb.workflows` submodule contains complete analysis workflows ...
+   For example, :mod:`alchemlyb.workflows.abfe` implements a complete absolute binding free energy calculation.".
 
 All of these components lend themselves well to writing clear and flexible pipelines for processing data needed for alchemical free energy calculations, and furthermore allow for scaling up via libraries like `dask`_ or `joblib`_.
 

docs/workflows.rst

@@ -7,10 +7,16 @@ of the results and step-by-step version that allows more flexibility.
 For developers, the skeleton of the workflow should follow the example in
 :class:`alchemlyb.workflows.base.WorkflowBase`.
 
+For users, **alchemlyb** offers a workflow :class:`alchemlyb.workflows.ABFE`
+similar to
+`Alchemical Analysis <https://github.com/MobleyLab/alchemical-analysis>`_
+for doing automatic absolute binding free energy (ABFE) analysis.
+
 .. currentmodule:: alchemlyb.workflows
 
 .. autosummary::
     :toctree: workflows
 
     base
+    ABFE
 

docs/workflows/alchemlyb.workflows.ABFE.rst

@@ -0,0 +1,161 @@
+The ABFE workflow
+==================
+The *Absolute binding free energy* (ABFE) workflow provides a complete workflow
+that uses the energy files generated by MD engine as input and generates the
+binding free energy as well as the analysis plots.
+
+Fully Automatic analysis
+------------------------
+*Absolute binding free energy* (ABFE) calculations can be analyzed with
+two lines of code in a fully automated manner (similar to 
+`Alchemical Analysis <https://github.com/MobleyLab/alchemical-analysis>`_).
+In this case, any parameters are set when invoking :class:`~alchemlyb.workflows.abfe.ABFE`
+and reasonable defaults are chosen for any parameters not set explicitly. The two steps
+are to
+
+1. initialize an instance of the :class:`~alchemlyb.workflows.abfe.ABFE` class
+2. invoke the :meth:`~alchemlyb.workflows.ABFE.run` method to execute
+   complete workflow.
+
+For a GROMACS ABFE simulation, executing the workflow would look similar
+to the following code (The log is configured by logger). ::
+
+    >>> from alchemtest.gmx import load_ABFE
+    >>> from alchemlyb.workflows import ABFE
+    >>> # Enable the logger
+    >>> import logging
+    >>> logging.basicConfig(filename='ABFE.log', level=logging.INFO)
+    >>> # Obtain the path of the data
+    >>> import os
+    >>> dir = os.path.dirname(load_ABFE()['data']['complex'][0])
+    >>> print(dir)
+    'alchemtest/gmx/ABFE/complex'
+    >>> workflow = ABFE(units='kcal/mol', software='Gromacs', dir=dir,
+    >>>                 prefix='dhdl', suffix='xvg', T=298, outdirectory='./')
+    >>> workflow.run(skiptime=10, uncorr='dhdl', threshold=50,
+    >>>              methods=('MBAR', 'BAR', 'TI'), overlap='O_MBAR.pdf',
+    >>>              breakdown=True, forwrev=10)
+
+
+The workflow uses the :class:`~alchemlyb.parsing` to parse the data from the
+energy files, remove the initial unequilibrated frames and decorrelate the data
+with :class:`~alchemlyb.preprocessing.subsampling`. The decorrelated dataset
+:ref:`dHdl <dHdl>` and :ref:`u_nk <u_nk>` are then passed to
+:class:`~alchemlyb.estimators` for free energy estimation. The workflow will
+also perform a set of analysis that allows the user to examine the quality of
+the estimation.
+
+File Input
+^^^^^^^^^^
+
+This command expects the energy files to be structured in two common ways. It
+could either be ::
+
+    simulation
+    ├── lambda_0
+    │   ├── prod.xvg
+    │   └── ...
+    ├── lambda_1
+    │   ├── prod.xvg
+    │   └── ...
+    └── ...
+
+Where :code:`dir='simulation/lambda_*', prefix='prod', suffix='xvg'`. Or ::
+
+    dhdl_files
+    ├── dhdl_0.xvg
+    ├── dhdl_1.xvg
+    └── ...
+
+Where :code:`dir='dhdl_files', prefix='dhdl_', suffix='xvg'`.
+
+Output
+^^^^^^
+
+The workflow returns the free energy estimate using all of
+:class:`~alchemlyb.estimators.TI`, :class:`~alchemlyb.estimators.BAR`,
+:class:`~alchemlyb.estimators.MBAR`. For ABFE calculations, the alchemical
+transformation is usually done is three stages, the *bonded*, *coul* and *vdw*
+which corresponds to the free energy contribution from applying the
+restraint to restrain the ligand to the protein, decouple/annihilate the
+coulombic interaction between the ligand and the protein and
+decouple/annihilate the protein-ligand lennard jones interactions. The result
+will be stored in :attr:`~alchemlyb.workflows.ABFE.summary` as
+:class:`pandas.Dataframe`. ::
+
+
+                          MBAR  MBAR_Error        BAR  BAR_Error         TI  TI_Error
+    States 0 -- 1     0.065967    0.001293   0.066544   0.001661   0.066663  0.001675
+           1 -- 2     0.089774    0.001398   0.089303   0.002101   0.089566  0.002144
+           2 -- 3     0.132036    0.001638   0.132687   0.002990   0.133292  0.003055
+    ...
+           26 -- 27   1.243745    0.011239   1.245873   0.015711   1.248959  0.015762
+           27 -- 28   1.128429    0.012859   1.124554   0.016999   1.121892  0.016962
+           28 -- 29   1.010313    0.016442   1.005444   0.017692   1.019747  0.017257
+    Stages coul      10.215658    0.033903  10.017838   0.041839  10.017854  0.048744
+           vdw       22.547489    0.098699  22.501150   0.060092  22.542936  0.106723
+           bonded     2.374144    0.014995   2.341631   0.005507   2.363828  0.021078
+           TOTAL     35.137291    0.103580  34.860619   0.087022  34.924618  0.119206
+
+Output Files
+^^^^^^^^^^^^
+
+For quality assessment, a couple of plots were generated and written to
+the folder specified by `outdirectory`.
+
+The :ref:`overlay matrix for the MBAR estimator <plot_overlap_matrix>` will be
+plotted and saved to :file:`O_MBAR.pdf`, which examines the overlap between
+different lambda windows.
+
+The :ref:`dHdl for TI <plot_TI_dhdl>` will be plotted to
+:file:`dhdl_TI.pdf`, allows one to examine if the lambda scheduling has
+covered the change of the gradient in the lambda space.
+
+The :ref:`dF states <plot_dF_states>` will be plotted to :file:`dF_state.pdf` in
+portrait model and :file:`dF_state_long.pdf` in landscape model, which
+allows the user to example the contributions from each lambda window.
+
+The forward and backward convergence will be plotted to :file:`dF_t.pdf` using
+:class:`~alchemlyb.estimators.MBAR` and saved in
+:attr:`~alchemlyb.workflows.ABFE.convergence`, which allows the user to
+examine if the simulation time is enough to achieve a converged result.
+
+Semi-automatic analysis
+-----------------------
+The same analysis could also performed in steps allowing access and modification
+to the data generated at each stage of the analysis. ::
+
+    >>> from alchemtest.gmx import load_ABFE
+    >>> from alchemlyb.workflows import ABFE
+    >>> # Obtain the path of the data
+    >>> import os
+    >>> dir = os.path.dirname(load_ABFE()['data']['complex'][0])
+    >>> print(dir)
+    'alchemtest/gmx/ABFE/complex'
+    >>> # Load the data
+    >>> workflow = ABFE(software='Gromacs', dir=dir,
+    >>>                 prefix='dhdl', suffix='xvg', T=298, outdirectory='./')
+    >>> # Set the unit.
+    >>> workflow.update_units('kcal/mol')
+    >>> # Read the data
+    >>> workflow.read()
+    >>> # Decorrelate the data.
+    >>> workflow.preprocess(skiptime=10, uncorr='dhdl', threshold=50)
+    >>> # Run the estimator
+    >>> workflow.estimate(methods=('mbar', 'bar', 'ti'))
+    >>> # Retrieve the result
+    >>> summary = workflow.generate_result()
+    >>> # Plot the overlap matrix
+    >>> workflow.plot_overlap_matrix(overlap='O_MBAR.pdf')
+    >>> # Plot the dHdl for TI
+    >>> workflow.plot_ti_dhdl(dhdl_TI='dhdl_TI.pdf')
+    >>> # Plot the dF states
+    >>> workflow.plot_dF_state(dF_state='dF_state.pdf')
+    >>> # Convergence analysis
+    >>> workflow.check_convergence(10, dF_t='dF_t.pdf')
+
+API Reference
+-------------
+.. autoclass:: alchemlyb.workflows.ABFE
+    :members:
+    :inherited-members:

src/alchemlyb/convergence/convergence.py

@@ -1,12 +1,15 @@
-import pandas as pd
 import logging
+from warnings import warn
+
+import pandas as pd
 import numpy as np
 
-from ..estimators import MBAR, BAR, TI, AutoMBAR
+from ..estimators import BAR, TI, FEP_ESTIMATORS, TI_ESTIMATORS
+from ..estimators import AutoMBAR as MBAR
 from .. import concat
 
 
-def forward_backward_convergence(df_list, estimator='mbar', num=10):
+def forward_backward_convergence(df_list, estimator='MBAR', num=10, **kwargs):
     '''Forward and backward convergence of the free energy estimate.
 
     Generate the free energy estimate as a function of time in both directions,
@@ -20,10 +23,16 @@ def forward_backward_convergence(df_list, estimator='mbar', num=10):
     ----------
     df_list : list
         List of DataFrame of either dHdl or u_nk.
-    estimator : {'mbar', 'bar', 'ti', 'autombar'}
+    estimator : {'MBAR', 'BAR', 'TI'}
         Name of the estimators.
+        See the important note below on the use of "MBAR".
+
+        .. deprecated:: 1.0.0
+           Lower case input is also accepted until release 2.0.0.
     num : int
         The number of time points.
+    kwargs : dict
+        Keyword arguments to be passed to the estimator.
 
     Returns
     -------
@@ -43,29 +52,40 @@ def forward_backward_convergence(df_list, estimator='mbar', num=10):
             9  3.044149       0.016405  3.044385        0.016402            1.0
 
 
+    Note
+    -----
+    :class:`~alchemlyb.estimators.AutoMBAR` is used when ``estimator='MBAR'``,
+    supply ``method`` keyword to restore the behavior of
+    :class:`~alchemlyb.estimators.MBAR`.
+    (:code:`forward_backward_convergence(u_nk, 'MBAR', num=2, method='adaptive')`)
+
+
     .. versionadded:: 0.6.0
+    .. versionchanged:: 1.0.0
+       The ``estimator`` accepts uppercase input.
+       The default for using ``estimator='MBAR'`` was changed from 
+       :class:`~alchemlyb.estimators.MBAR` to :class:`~alchemlyb.estimators.AutoMBAR`.
 
     '''
     logger = logging.getLogger('alchemlyb.convergence.'
                                'forward_backward_convergence')
     logger.info('Start convergence analysis.')
     logger.info('Check data availability.')
+    if estimator.upper() != estimator:
+        warn("Using lower-case strings for the 'estimator' kwarg in "
+             "convergence.forward_backward_convergence() is deprecated in "
+             "1.0.0 and only upper case will be accepted in 2.0.0",
+            DeprecationWarning)
+        estimator = estimator.upper()
 
-    if estimator.lower() == 'mbar':
-        logger.info('Use MBAR estimator for convergence analysis.')
-        estimator_fit = MBAR().fit
-    elif estimator.lower() == 'autombar':
-        logger.info('Use AutoMBAR estimator for convergence analysis.')
-        estimator_fit = AutoMBAR().fit
-    elif estimator.lower() == 'bar':
-        logger.info('Use BAR estimator for convergence analysis.')
-        estimator_fit = BAR().fit
-    elif estimator.lower() == 'ti':
-        logger.info('Use TI estimator for convergence analysis.')
-        estimator_fit = TI().fit
+    if estimator not in (FEP_ESTIMATORS + TI_ESTIMATORS):
+        msg = f"Estimator {estimator} is not available in {FEP_ESTIMATORS + TI_ESTIMATORS}."
+        logger.error(msg)
+        raise ValueError(msg)
     else:
-        raise ValueError(
-            '{} is not a valid estimator.'.format(estimator))
+        # select estimator class by name
+        estimator_fit = globals()[estimator](**kwargs).fit
+        logger.info(f'Use {estimator} estimator for convergence analysis.')
 
     logger.info('Begin forward analysis')
     forward_list = []

src/alchemlyb/estimators/__init__.py

@@ -1,3 +1,6 @@
 from .mbar_ import MBAR, AutoMBAR
 from .bar_ import BAR
 from .ti_ import TI
+
+FEP_ESTIMATORS = [MBAR.__name__, AutoMBAR.__name__, BAR.__name__]
+TI_ESTIMATORS = [TI.__name__]

src/alchemlyb/estimators/mbar_.py

@@ -152,6 +152,17 @@ class AutoMBAR(MBAR):
     :class:`AutoMBAR` may be useful in high-throughput calculations where it can avoid 
     failures due non-converged MBAR estimates.
 
+    Parameters
+    ----------
+
+    method : str, optional, default=None
+        The optimization routine to use. This parameter defaults to ``None``.
+        When a specific method is set, AutoMBAR will behave in the same way
+        as MBAR.
+
+        .. versionadded:: 1.0.0
+
+
     Note
     ----
     All arguments are described under :class:`MBAR` except that the solver method 
@@ -163,30 +174,35 @@ class AutoMBAR(MBAR):
 
 
     .. versionadded:: 0.6.0
+    .. versionchanged:: 1.0.0
+       AutoMBAR accepts the `method` argument.
     """
     def __init__(self, maximum_iterations=10000, relative_tolerance=1.0e-7,
-                 initial_f_k=None, verbose=False):
+                 initial_f_k=None, verbose=False, method=None):
         super().__init__(maximum_iterations=maximum_iterations,
                          relative_tolerance=relative_tolerance,
                          initial_f_k=initial_f_k,
-                         verbose=verbose, method=None)
+                         verbose=verbose, method=method)
         self.logger = logging.getLogger('alchemlyb.estimators.AutoMBAR')
 
     def _do_MBAR(self, u_nk, N_k, solver_protocol):
-        self.logger.info('Initialise the automatic routine of the MBAR '
-                         'estimator.')
-        # Try the fastest method first
-        try:
-            self.logger.info('Trying the hybr method.')
-            solver_protocol["method"] = 'hybr'
-            mbar, out = super()._do_MBAR(u_nk, N_k, solver_protocol)
-        except pymbar.utils.ParameterError:
+        if solver_protocol["method"] is None:
+            self.logger.info('Initialise the automatic routine of the MBAR '
+                             'estimator.')
+            # Try the fastest method first
             try:
-                self.logger.info('Trying the adaptive method.')
-                solver_protocol["method"] = 'adaptive'
+                self.logger.info('Trying the hybr method.')
+                solver_protocol["method"] = 'hybr'
                 mbar, out = super()._do_MBAR(u_nk, N_k, solver_protocol)
             except pymbar.utils.ParameterError:
-                self.logger.info('Trying the BFGS method.')
-                solver_protocol["method"] = 'BFGS'
-                mbar, out = super()._do_MBAR(u_nk, N_k, solver_protocol)
-        return mbar, out
+                try:
+                    self.logger.info('Trying the adaptive method.')
+                    solver_protocol["method"] = 'adaptive'
+                    mbar, out = super()._do_MBAR(u_nk, N_k, solver_protocol)
+                except pymbar.utils.ParameterError:
+                    self.logger.info('Trying the BFGS method.')
+                    solver_protocol["method"] = 'BFGS'
+                    mbar, out = super()._do_MBAR(u_nk, N_k, solver_protocol)
+            return mbar, out
+        else:
+            return super()._do_MBAR(u_nk, N_k, solver_protocol)