diff --git a/package/AUTHORS b/package/AUTHORS index ed82321a9e..3ed92cc27b 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -284,6 +284,7 @@ Chronological list of authors - Sai Udayagiri - Apoorva Verma - Aryaman Chaudhri + - Srikar JY External code ------------- diff --git a/package/CHANGELOG b/package/CHANGELOG index c043f683f0..c9c5e34b11 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -18,11 +18,12 @@ The rules for this file: spyke7, talagayev, tanii1125, BradyAJohnston, hejamu, jeremyleung521, harshitgajjela-droid, kunjsinha, aygarwal, jauy123, Dreamstick9, ollyfutur, Amarendra22, charity-g, ParthUppal523, apoorva-01, RMeli, - raulloiscuns, Aryaman-Chaudhri + raulloiscuns, Aryaman-Chaudhri, srikarjy * 2.11.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Added `.gitattributes` to enforce LF (\n) as line endings and renormalized existing files to conform (Issue #5315, PR #5446) * `AtomGroup.rotate()` and the `rotateby` trajectory transformation now diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 7a0510eab0..70fea80fa0 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -316,8 +316,13 @@ def fit(self): return self.results.fit def _conclude(self): + if self.n_frames == 0: + raise ValueError( + "PersistenceLength.run() analyzed zero frames; check the " + "start/stop/step or frames selection passed to run()." + ) norm = np.linspace(self.chainlength - 1, 1, self.chainlength - 1) - norm *= len(self._atomgroups) * self._trajectory.n_frames + norm *= len(self._atomgroups) * self.n_frames self.results.bond_autocorrelation = ( self.results.raw_bond_autocorr / norm ) diff --git a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py index 54005c8a66..360d4c6651 100644 --- a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py +++ b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py @@ -31,7 +31,7 @@ import matplotlib import matplotlib.pyplot as plt -from numpy.testing import assert_almost_equal, assert_equal +from numpy.testing import assert_allclose, assert_almost_equal, assert_equal from MDAnalysisTests.datafiles import Plength, TRZ_psf, TRZ @@ -177,3 +177,218 @@ def test_circular(self): u.add_TopologyAttr(Bonds(bondlist)) with pytest.raises(ValueError, match="Cyclical"): polymer.sort_backbone(u.atoms) + + +def _reference_bond_autocorrelation(chains, frame_indices): + """Independently compute the expected normalized bond autocorrelation + for exactly the given sequence of trajectory frame indices. + + Mirrors the mathematical definition directly from atom positions, + without calling into :class:`~MDAnalysis.analysis.polymer. + PersistenceLength` or reusing any of its results, so this stays a + meaningful check even if the class under test is itself broken. + Frame indices are visited explicitly (in the order given, duplicates + allowed) rather than via a slice, so this can independently mirror + any of `run`'s `start`/`stop`/`step`/`frames` selections. See Issue + #5453. + """ + chainlength = len(chains[0]) + raw = np.zeros(chainlength - 1, dtype=np.float64) + + universe = chains[0].universe + n_frames = 0 + for frame in frame_indices: + universe.trajectory[frame] + n_frames += 1 + for chain in chains: + vecs = chain.positions[1:] - chain.positions[:-1] + vecs = vecs / np.sqrt((vecs * vecs).sum(axis=1))[:, None] + inner_pr = np.inner(vecs, vecs) + for i in range(chainlength - 1): + raw[: (chainlength - 1) - i] += inner_pr[i, i:] + + norm = np.linspace(chainlength - 1, 1, chainlength - 1) + norm *= len(chains) * n_frames + return raw / norm, n_frames + + +class TestPersistenceLengthSlicedNormalization(object): + # Regression test for Issue #5453. + @staticmethod + @pytest.fixture(scope="class") + def u(): + return mda.Universe(TRZ_psf, TRZ) + + @staticmethod + @pytest.fixture(scope="class") + def chains(u): + backbones = [ + chain.select_atoms("not name O* H*") for chain in u.atoms.fragments + ] + return [polymer.sort_backbone(bb) for bb in backbones] + + def test_full_trajectory_unchanged(self, chains, u): + p = polymer.PersistenceLength(chains).run() + + assert p.n_frames == u.trajectory.n_frames + + expected, n_frames = _reference_bond_autocorrelation( + chains, range(u.trajectory.n_frames) + ) + assert n_frames == u.trajectory.n_frames + # atol/rtol (looser than a plain decimal=6 comparison) account for + # PersistenceLength accumulating raw_bond_autocorr in float32 while + # this reference accumulates in float64; their rounding drift is + # a couple of float32 ULPs at this trajectory's frame/chain count, + # not a sign of disagreement. + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + ) + + def test_sliced_run_normalized_by_frames_analyzed(self, chains, u): + n_sliced_frames = 3 + assert n_sliced_frames < u.trajectory.n_frames + + p = polymer.PersistenceLength(chains).run(stop=n_sliced_frames) + + assert p.n_frames == n_sliced_frames + assert p.n_frames != u.trajectory.n_frames + + expected, n_frames = _reference_bond_autocorrelation( + chains, range(n_sliced_frames) + ) + assert n_frames == n_sliced_frames + # atol/rtol (looser than a plain decimal=6 comparison) account for + # PersistenceLength accumulating raw_bond_autocorr in float32 while + # this reference accumulates in float64; their rounding drift is + # a couple of float32 ULPs at this trajectory's frame/chain count, + # not a sign of disagreement. + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + ) + + # Sanity check: reproduce the old (buggy) behavior by normalizing + # the *sliced* raw sum with the *total* trajectory frame count + # instead of n_sliced_frames, and confirm that no longer matches + # the class's output -- i.e. this test is sensitive to the bug + # described in Issue #5453. + wrong = expected * n_sliced_frames / u.trajectory.n_frames + assert not np.allclose(p.results.bond_autocorrelation, wrong) + + @pytest.mark.parametrize( + "run_kwargs, frame_indices", + [ + pytest.param({"step": 2}, [0, 2, 4], id="step_only"), + pytest.param( + {"start": 1, "stop": 6, "step": 2}, + [1, 3, 5], + id="start_stop_step", + ), + pytest.param({"frames": [0, 2, 5]}, [0, 2, 5], id="frames_list"), + pytest.param( + {"frames": [4, 1, 3]}, + [4, 1, 3], + id="frames_list_unordered", + ), + ], + ) + def test_various_frame_selections_normalized_correctly( + self, chains, u, run_kwargs, frame_indices + ): + # Regression test for Issue #5453: the normalization fix relies + # on self.n_frames, which AnalysisBase derives the same way + # regardless of whether the run was sliced via start/stop/step + # or via an explicit frames= list. Exercise more than just a + # plain stop= slice to guard against that assumption breaking. + p = polymer.PersistenceLength(chains).run(**run_kwargs) + + assert p.n_frames == len(frame_indices) + assert p.n_frames != u.trajectory.n_frames + + expected, n_frames = _reference_bond_autocorrelation( + chains, frame_indices + ) + assert n_frames == len(frame_indices) + # atol/rtol (looser than a plain decimal=6 comparison) account for + # PersistenceLength accumulating raw_bond_autocorr in float32 while + # this reference accumulates in float64; their rounding drift is + # a couple of float32 ULPs at this trajectory's frame/chain count, + # not a sign of disagreement. + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + ) + + def test_sliced_run_parallel_matches_serial(self, chains, u): + # Regression test for Issue #5453: the bug was introduced by the + # multi-worker refactor in PR #5074, but neither existing test + # actually exercised backend="multiprocessing" together with a + # sliced run. _conclude() only ever runs once, on the main + # process, after worker results are merged -- so self.n_frames + # at normalization time should be the total sliced frame count + # regardless of how work was chunked across workers. Confirm + # that holds, rather than only checking it by hand. + n_sliced_frames = 5 + assert n_sliced_frames < u.trajectory.n_frames + + serial = polymer.PersistenceLength(chains).run(stop=n_sliced_frames) + parallel = polymer.PersistenceLength(chains).run( + stop=n_sliced_frames, backend="multiprocessing", n_workers=2 + ) + + assert serial.n_frames == parallel.n_frames == n_sliced_frames + assert_allclose( + serial.results.bond_autocorrelation, + parallel.results.bond_autocorrelation, + rtol=1e-4, + atol=5e-6, + ) + + def test_empty_frame_selection_raises_clear_error(self, chains): + # Edge case found while verifying the fix for Issue #5453: an + # empty frame selection drives self.n_frames to 0, which used + # to silently divide by the (nonzero) total trajectory frame + # count and return a fabricated all-zero curve with no warning. + # After the fix, the zero-frame divisor would instead produce + # NaNs that surface as an opaque scipy error deep inside + # curve_fit. Raise a clear, immediate error instead. + with pytest.raises(ValueError, match="analyzed zero frames"): + polymer.PersistenceLength(chains).run(start=2, stop=2) + + def test_single_atom_chain_unaffected_by_normalization_fix(self, u): + # Not a regression from Issue #5453's fix: a chain with only one + # atom has no bonds to correlate (chainlength - 1 == 0), so + # raw_bond_autocorr and norm are both length-0 arrays regardless + # of the run() call's start/stop/step. This has always raised + # from inside _perform_fit()'s curve_fit call, independent of + # which frame count _conclude() normalizes by. Documented here + # so it doesn't get mistaken for something this PR should fix. + one_atom_chain = [u.atoms.fragments[0][:1]] + with pytest.raises(ValueError): + polymer.PersistenceLength(one_atom_chain).run() + + @pytest.mark.xfail( + reason=( + "Pre-existing AnalysisBase bug, independent of Issue #5453: " + "_setup_computation_groups checks `isinstance(obj, bool)`, " + "which is False for numpy bool scalars, so a numpy boolean " + "frames= mask is treated as an integer index array (True/" + "False read as 1/0) instead of a mask. A plain Python list " + "of bools is unaffected. Left xfail to document the " + "framework-level issue without attempting to fix " + "AnalysisBase in this PR." + ), + strict=True, + ) + def test_numpy_bool_mask_frames_selection(self, chains, u): + mask = np.zeros(u.trajectory.n_frames, dtype=bool) + mask[[0, 2, 4]] = True + + p = polymer.PersistenceLength(chains).run(frames=mask) + + assert p.n_frames == int(mask.sum()) + expected, n_frames = _reference_bond_autocorrelation( + chains, np.flatnonzero(mask) + ) + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + )