Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

replace scipy.integrate.simps with simpson #282

Merged
merged 6 commits into from
Oct 10, 2024
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 3 additions & 6 deletions .github/workflows/ci.yaml
Original file line number Diff line number Diff line change
@@ -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 }}
9 changes: 6 additions & 3 deletions CHANGES
Original file line number Diff line number Diff line change
@@ -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()
4 changes: 2 additions & 2 deletions devtools/conda-envs/test_env.yaml
Original file line number Diff line number Diff line change
@@ -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
5 changes: 2 additions & 3 deletions mdpow/equil.py
Original file line number Diff line number Diff line change
@@ -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...
18 changes: 14 additions & 4 deletions mdpow/fep.py
Original file line number Diff line number Diff line change
@@ -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[
10 changes: 8 additions & 2 deletions mdpow/tests/test_analysis.py
Original file line number Diff line number Diff line change
@@ -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):
3 changes: 1 addition & 2 deletions setup.py
Original file line number Diff line number Diff line change
@@ -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",