replace scipy.integrate.simps with simpson

.github/workflows/ci.yaml

@@ -34,9 +34,6 @@ jobs:
         python-version: ["3.10"]
         gromacs-version: ["4.6.5", "2018.6", "2020.6", "2021.1", "2022.4", "2023.1"]
         include:
-          - os: ubuntu-latest
-            python-version: "3.8"
-            gromacs-version: "2023.1"
           - os: ubuntu-latest
             python-version: "3.9"
             gromacs-version: "2023.1"
@@ -49,10 +46,10 @@ jobs:
 
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: mamba environment and package installation
-      uses: mamba-org/setup-micromamba@v1
+      uses: mamba-org/setup-micromamba@v2
       with:
         environment-file: devtools/conda-envs/test_env.yaml
         condarc: |
@@ -93,7 +90,7 @@ jobs:
         pytest -v --durations=20 --cov=mdpow --cov-report=xml --color=yes ./mdpow/tests
 
     - name: Codecov
-      uses: codecov/codecov-action@v3
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         name: codecov-${{ matrix.os }}-py${{ matrix.python-version }}

CHANGES

@@ -5,14 +5,17 @@ CHANGES for MDPOW
 Add summary of changes for each release. Use ISO 8061 dates. Reference
 GitHub issues numbers and PR numbers.
 
-2023-??-??    0.9.0
+2024-??-??    0.9.0
 cadeduckworth, orbeckst, VOD555, a-ws-m
 
 Changes
 
+* change in TI implementation in fep.Gsolv.analysis(): scipy.integrate.simpson()
+  now always uses Cartwright's approach to compute the last interval instead of
+  the old `even="last"` behavior. This change **may lead to small numerical 
+  differences in output** (#281)
 * added support for Python 3.10 (#202)
-* dropped testing on Python 3.6 (PR #220, #202)
-* dropped testing on Python 3.7 (minimally supported Python >= 3.8, #248)
+* dropped testing/support for Python 3.8 (#281), 3.7 (#248). 3.6 (PR #220, #202)
 * support Gromacs 2022.4 and 2023.1 (#256)
 * use pymbar >= 4 and alchemlyb >= 2 (#246)
 * for ensemble.EnsembleAnalysis._single_frame() 

devtools/conda-envs/test_env.yaml

@@ -6,7 +6,7 @@ dependencies:
 - python
 - six
 - numpy
-- scipy
+- scipy >=1.11.0
 - matplotlib-base
 - pandas
 - scikit-learn
@@ -22,7 +22,7 @@ dependencies:
 - cairosvg
 - pypdf
 
-  # Testing
+# Testing
 - pytest
 - pytest-pep8
 - pytest-cov

mdpow/equil.py

@@ -535,9 +535,8 @@ def _MD(self, protocol, **kwargs):
         kwargs["top"] = self.files.topology
         kwargs["includes"] = asiterable(kwargs.pop("includes", [])) + self.dirs.includes
         kwargs["ndx"] = self.files.ndx
-        kwargs[
-            "mainselection"
-        ] = None  # important for SD (use custom mdp and ndx!, gromacs.setup._MD)
+        # important for SD (use custom mdp and ndx!, gromacs.setup._MD):
+        kwargs["mainselection"] = None
         self._checknotempty(kwargs["struct"], "struct")
         if not os.path.exists(kwargs["struct"]):
             # struct is not reliable as it depends on qscript so now we just try everything...

mdpow/fep.py

@@ -1020,7 +1020,7 @@ def analyze(self, force=False, stride=None, autosave=True, ncorrel=25000):
 
         The dV/dlambda graphs are integrated with the composite Simpson's rule
         (and if the number of datapoints are even, the first interval is
-        evaluated with the trapezoidal rule); see :func:`scipy.integrate.simps`
+        evaluated with the trapezoidal rule); see :func:`scipy.integrate.simpson`
         for details). Note that this implementation of Simpson's rule does not
         require equidistant spacing on the lambda axis.
 
@@ -1081,6 +1081,14 @@ def analyze(self, force=False, stride=None, autosave=True, ncorrel=25000):
                              Diego 2002
 
         .. _p526: http://books.google.co.uk/books?id=XmyO2oRUg0cC&pg=PA526
+
+        .. versionchanged:: 0.9.0
+           Change in how the Simpson's rule integration algorithm (namely,
+           :func:`scipy.integrate.simpson`) handles even number of intervals:
+           Previously, the old `even="last"` was used but since scipy 1.11.0
+           Cartwright's approach is used. This change **leads to numerically
+           slightly different results** between MDPOW 0.9.0 and earlier
+           versions.
         """
         stride = stride or self.stride
         logger.info("Analysis stride is %s.", stride)
@@ -1137,9 +1145,11 @@ def analyze(self, force=False, stride=None, autosave=True, ncorrel=25000):
                 "tcorrel": tc,
             }
             # Combined Simpson rule integration:
-            # even="last" because dV/dl is smoother at the beginning so using trapezoidal
-            # integration there makes less of an error (one hopes...)
-            a = scipy.integrate.simps(Y, x=lambdas, even="last")
+            # Used to have 'even="last"' because dV/dl is smoother at the beginning so
+            # using trapezoidal integration there makes less of an error (one hopes...)
+            # but recent versions of scipy (eg 1.14) always use Cartwright's approach
+            # for the last interval and "even" is not a kwarg anymore.
+            a = scipy.integrate.simpson(Y, x=lambdas)
             da = numkit.integration.simps_error(DY, x=lambdas, even="last")
             self.results.DeltaA[component] = QuantityWithError(a, da)
             GibbsFreeEnergy += self.results.DeltaA[

mdpow/tests/test_analysis.py

@@ -84,17 +84,23 @@ def assert_DeltaA(G):
         #   original values are only reproduced to 5 decimals, see PR #166"
         # - June 2023: in CI, >= 3.8 results differ from reference values (although
         #   locally no changes are obvious) after ~4 decimals for unknown reasons.
+        # - Oct 2024: change to scipy.integrate.simpson(): use Cartwright's approach
+        #   instead of even="last": changes the mean (DeltaA: from -3.722 to now -3.643)
         DeltaA = G.results.DeltaA
         assert_array_almost_equal(
-            DeltaA.Gibbs.astuple(), (-3.7217472974883794, 2.3144288928034911), decimal=3
+            DeltaA.Gibbs.astuple(),
+            (-3.6429995060434432, 2.3141470255028795),
+            decimal=3,
         )
         assert_array_almost_equal(
             DeltaA.coulomb.astuple(),
             (8.3346255170099575, 0.73620918517131495),
             decimal=3,
         )
         assert_array_almost_equal(
-            DeltaA.vdw.astuple(), (-4.6128782195215781, 2.1942144688960972), decimal=3
+            DeltaA.vdw.astuple(),
+            (-4.691626010966514, 2.1940230979105584),
+            decimal=3,
         )
 
     def test_convert_edr(self, fep_benzene_directory):

setup.py

@@ -27,7 +27,6 @@
         "Operating System :: MacOS :: MacOS X",
         "Programming Language :: Python",
         "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.8",
         "Programming Language :: Python :: 3.9",
         "Programming Language :: Python :: 3.10",
         "Topic :: Scientific/Engineering :: Chemistry",
@@ -60,7 +59,7 @@
     },
     install_requires=[
         "numpy>=1.6",
-        "scipy",
+        "scipy>=1.11.0",
         "pyyaml",
         "GromacsWrapper>=0.5.1",
         "numkit",