Update statistical inefficiency calculation with statistical_inefficiency in alchemlyb 2.1.0

.github/workflows/ci.yaml

@@ -77,8 +77,10 @@ jobs:
       run: |
          python -m pip install pybol
     - name: Install GROMACS (${{ matrix.gromacs-version }})
+      # include mdanalysis here AGAIN to work around micromamba removing mdanalysis 
+      # while downgrading libxml2 (see issue #252)
       run: |
-         micromamba install 'gromacs==${{ matrix.gromacs-version }}' pocl
+         micromamba install 'gromacs==${{ matrix.gromacs-version }}' pocl mdanalysis
 
     - name: Install package (with no dependencies)
       run: |

CHANGES

@@ -18,9 +18,10 @@ Changes
   and use _single_frame OR _single_universe (#216)
 * _prepare_universe and _conclude_universe removed from 
   EnsembleAnalysis.run() method, no longer needed (per comments, #199)
-* added support for Python 3.10
-* dropped testing on Python 3.6
-* internal log_banner() now uses logger as argument
+* added support for Python 3.10 (#202)
+* dropped testing on Python 3.6 (PR #220, #202)
+* use pymbar >= 4 and alchemlyb >= 2 (#246)
+* internal log_banner() now uses logger as argument (PR #247)
 
 Enhancements
 

devtools/conda-envs/test_env.yaml

@@ -14,7 +14,8 @@ dependencies:
 - mdanalysis >=2
 - numkit
 - gromacswrapper
-- alchemlyb <2
+- alchemlyb >=2
+- pymbar >=4
 - rdkit
 - seaborn
 

doc/requirements.txt

@@ -4,7 +4,8 @@ scipy
 pandas
 pyyaml
 GromacsWrapper>=0.5.1
-alchemlyb
+alchemlyb>=2
+pymbar>=4
 mdanalysis
 rdkit
 seaborn

mdpow/fep.py

@@ -157,11 +157,10 @@
 from numkit.observables import QuantityWithError
 
 from alchemlyb.parsing.gmx import extract_dHdl, extract_u_nk
-from alchemlyb.estimators import TI, BAR
-from alchemlyb.estimators import AutoMBAR as MBAR
+from alchemlyb.estimators import TI, BAR, MBAR
 from alchemlyb.parsing.gmx import _extract_dataframe
-from pymbar.timeseries import (statisticalInefficiency,
-                               subsampleCorrelatedData, )
+from alchemlyb.preprocessing.subsampling import statistical_inefficiency
+
 import gromacs
 import gromacs.utilities
 try:
@@ -1054,19 +1053,16 @@ def collect_alchemlyb(self, SI=True, start=0, stop=None, stride=None, autosave=T
             for l in lambdas:
                 xvg_file = self.dgdl_xvg(self.wdir(component, l))
                 xvg_df = extract(xvg_file, T=self.Temperature).iloc[start:stop:stride]
+                full_len = len(xvg_df)
                 if SI:
                     logger.info("Performing statistical inefficiency analysis for window %s %04d" % (component, 1000 * l))
                     ts = _extract_dataframe(xvg_file).iloc[start:stop:stride]
                     ts = pd.DataFrame({'time': ts.iloc[:,0], 'dhdl': ts.iloc[:,1]})
                     ts = ts.set_index('time')
-                    # calculate statistical inefficiency of series
-                    statinef  = statisticalInefficiency(ts, fast=False)
-                    logger.info("The statistical inefficiency value is {:.4f}.".format(statinef))
-                    logger.info("The data are subsampled every {:d} frames.".format(int(np.ceil(statinef))))
-                    # use the subsampleCorrelatedData function to get the subsample index
-                    indices = subsampleCorrelatedData(ts, g=statinef,
-                                                      conservative=True)
-                    xvg_df = xvg_df.iloc[indices]
+                    # use the statistical_inefficiency function to subsample the data
+                    xvg_df = statistical_inefficiency(xvg_df, ts, conservative=True)
+                    logger.info("The statistical inefficiency value is {:.4f}.".format(full_len/len(xvg_df)/2))
+                    logger.info("The data are subsampled every {:d} frames.".format(int(np.ceil(full_len/len(xvg_df)/2))))
                 val.append(xvg_df)
             self.results.xvg[component] = (np.array(lambdas), pd.concat(val))
 

setup.py

@@ -57,9 +57,9 @@
                         'numkit',
                         'six',
                         'mdanalysis>=2',
-                        'alchemlyb<2',
+                        'alchemlyb>=2',
                         'pandas',
-                        'pymbar',
+                        'pymbar>=4',
                         'matplotlib',
                         'seaborn',
                         'rdkit',