replace scipy.integrate.simps with simpson

Becksteinlab · Oct 10, 2024 · a5162d8 · a5162d8
1 parent 120c5dd
commit a5162d8
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 7 deletions.
diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml
@@ -6,7 +6,7 @@ dependencies:
 - python
 - six
 - numpy
-- scipy
+- scipy >=1.6.0
 - matplotlib-base
 - pandas
 - scikit-learn
@@ -22,7 +22,7 @@ dependencies:
 - cairosvg
 - pypdf
 
-  # Testing
+# Testing
 - pytest
 - pytest-pep8
 - pytest-cov

diff --git a/mdpow/fep.py b/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.
 
@@ -1137,9 +1137,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[

diff --git a/setup.py b/setup.py
@@ -60,7 +60,7 @@
     },
     install_requires=[
         "numpy>=1.6",
-        "scipy",
+        "scipy>=1.6.0",
         "pyyaml",
         "GromacsWrapper>=0.5.1",
         "numkit",