Skip to content

Fix compatibility with numpy 2.5 #5404

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.

Sign up for GitHub

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

Merged
RMeli merged 7 commits into from
Jun 23, 2026
+54 −32
Merged

Fix compatibility with numpy 2.5 #5404

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0df2fd7
Fix compatibility with numpy 2.5
IAlibay Jun 22, 2026
a0dce8c
Merge branch 'develop' into fix-ci
IAlibay Jun 22, 2026
3bcf762
fix CI
IAlibay Jun 22, 2026
fd6fa86
format two files
IAlibay Jun 22, 2026
7cecd60
Try something a little bit more clever to compare the arrays
IAlibay Jun 22, 2026
0aebde9
add more comments
IAlibay Jun 22, 2026
2f46b14
Update package/CHANGELOG
orbeckst Jun 22, 2026
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
4 changes: 4 additions & 0 deletions package/CHANGELOG
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,10 @@ The rules for this file:
* 2.11.0

Fixes
* The `principal_axes` method in :class:`Masses` now uses
`np.linalg.eigh` instead of `np.linalg.eig`, improving numerical
stability. This may lead to slightly different results from previous
version of MDAnalysis (#5403, PR #5404).
* `MDAnalysis.analysis.nucleicacids.WatsonCrickDist`, `MinorPairDist`,
and `MajorPairDist` now match residue names against the full resname
instead of only the first character, fixing incorrect behaviour with
Expand Down
2 changes: 1 addition & 1 deletion package/MDAnalysis/core/topologyattrs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2084,7 +2084,7 @@ def principal_axes(group, wrap=False):
is deprecated and will be removed in version 3.0.
"""
atomgroup = group.atoms
e_val, e_vec = np.linalg.eig(atomgroup.moment_of_inertia(wrap=wrap))
e_val, e_vec = np.linalg.eigh(atomgroup.moment_of_inertia(wrap=wrap))

@IAlibay IAlibay Jun 22, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this should fix it?

Someone with better maths knowledge should double check this - my understanding is that the moment of intertia tensor is always a real symmetric matrix, so it should satisfy the requirements of np.linalg.eigh.

My understanding is that eigh doesn't have the stability issues of eig, so I tried swapping it out and it seems to be working fine.

@RMeli RMeli Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed, eigh is preferred for symmetric/Hermitian matrices.


# Sort
indices = np.argsort(e_val)[::-1]
Expand Down
39 changes: 26 additions & 13 deletions testsuite/MDAnalysisTests/core/test_atomgroup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1280,8 +1280,8 @@ class TestPBCFlag(object):
),
"principal_axes": np.array(
[
[0.78787867, 0.26771575, -0.55459488],
[-0.40611024, -0.45112859, -0.7947059],
[-0.78787867, -0.26771575, 0.55459488],

@IAlibay IAlibay Jun 22, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe these are fine since the eigenvector sign is arbitrary?

@RMeli RMeli Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the sign is algorithm-dependent and therefore can differ between eig and eigh.

@RMeli RMeli Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like some tests are still failing because of sign issues. Maybe to make them more robust we could check for all the eigenvector tests that they are numerically equivalent up to a sign change?

@IAlibay IAlibay Jun 22, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I'm thinking we should do an abs check, will update the tests accordingly.

@orbeckst orbeckst Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good idea

@IAlibay IAlibay Jun 22, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done - just to sanity check me, can I get a re-review please?

@orbeckst orbeckst Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought about it again. abs is not the same as sign-check. Will comment on review.

[0.40611024, 0.45112859, 0.7947059],
[-0.46294889, 0.85135849, -0.24671249],
]
),
Expand Down Expand Up @@ -1315,8 +1315,8 @@ class TestPBCFlag(object):
),
"principal_axes": np.array(
[
[0.85911708, -0.19258726, -0.4741603],
[0.07520116, 0.96394227, -0.25526473],
[-0.85911708, 0.19258726, 0.4741603],
[-0.07520116, -0.96394227, 0.25526473],
[0.50622389, 0.18364489, 0.84262206],
]
),
Expand Down Expand Up @@ -1355,6 +1355,14 @@ def test_wrap(self, ag, wrap, ref, method_name):
if method_name == "bsphere":
assert_almost_equal(result[0], ref[method_name][0], self.prec)
assert_almost_equal(result[1], ref[method_name][1], self.prec)
elif method_name == "principal_axes":
# See PR #5404
# The direction (sign) of the principal axes is dependent on the
# specific algorithm used, but the direction itself is not physically
# relevant, so we get the signs to flip any anti-parallel vectors before
# comparing the two results arrays
signs = np.sign(np.einsum("ij,ij->i", result, ref[method_name]))

@orbeckst orbeckst Jun 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fancy, I convinced myself that it does the right thing for the output of principal_axes().

assert_almost_equal(result * signs[:, np.newaxis], ref[method_name], self.prec)
else:
assert_almost_equal(result, ref[method_name], self.prec)

Expand Down Expand Up @@ -1620,16 +1628,21 @@ def test_coordinates(self, ag):
)

def test_principal_axes(self, ag):
assert_almost_equal(
ag.principal_axes(),
np.array(
[
[1.53389276e-03, 4.41386224e-02, 9.99024239e-01],
[1.20986911e-02, 9.98951474e-01, -4.41539838e-02],
[-9.99925632e-01, 1.21546132e-02, 9.98264877e-04],
]
),
ref = np.array(
[
[-1.53389276e-03, -4.41386224e-02, -9.99024239e-01],
[-1.20986911e-02, -9.98951474e-01, 4.41539838e-02],
[-9.99925632e-01, 1.21546132e-02, 9.98264877e-04],
]
)
result = ag.principal_axes()
# See PR #5404
# The direction (sign) of the principal axes is dependent on the
# specific algorithm used, but the direction itself is not physically
# relevant, so we get the signs to flip any anti-parallel vectors before
# comparing the two results arrays
signs = np.sign(np.einsum("ij,ij->i", result, ref))
assert_almost_equal(result * signs[:, np.newaxis], ref)

def test_principal_axes_duplicates(self, ag):
ag2 = ag + ag[0]
Expand Down
18 changes: 9 additions & 9 deletions testsuite/MDAnalysisTests/core/test_groups.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,16 +243,16 @@ def test_passive_decorator(self, ag):
assert_almost_equal(ag.radius_of_gyration(), 2.400527938286)
assert_almost_equal(ag.shape_parameter(), 0.61460819)
assert_almost_equal(ag.asphericity(), 0.4892751412)
assert_almost_equal(
ag.principal_axes(),
np.array(
[
[0.7574113, -0.113481, 0.643001],
[0.5896252, 0.5419056, -0.5988993],
[-0.2804821, 0.8327427, 0.4773566],
]
),
ref_pa = np.array(
[
[-0.7574113, 0.113481, -0.643001],
[-0.5896252, -0.5419056, 0.5988993],
[-0.2804821, 0.8327427, 0.4773566],
]
)
result_pa = ag.principal_axes()
signs_pa = np.sign(np.einsum("ij,ij->i", result_pa, ref_pa))
assert_almost_equal(result_pa * signs_pa[:, np.newaxis], ref_pa)
assert_almost_equal(
ag.center_of_charge(),
np.array([11.0800112, 8.8885659, -8.9886632]),
Expand Down
23 changes: 14 additions & 9 deletions testsuite/MDAnalysisTests/core/test_topologyattrs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -417,16 +417,21 @@ def ag(self):
return universe.atoms # prototypical AtomGroup

def test_principal_axes(self, ag):
assert_almost_equal(
ag.principal_axes(),
np.array(
[
[1.53389276e-03, 4.41386224e-02, 9.99024239e-01],
[1.20986911e-02, 9.98951474e-01, -4.41539838e-02],
[-9.99925632e-01, 1.21546132e-02, 9.98264877e-04],
]
),
ref = np.array(
[
[-1.53389276e-03, -4.41386224e-02, -9.99024239e-01],
[-1.20986911e-02, -9.98951474e-01, 4.41539838e-02],
[-9.99925632e-01, 1.21546132e-02, 9.98264877e-04],
]
)
result = ag.principal_axes()
# See PR #5404
# The direction (sign) of the principal axes is dependent on the
# specific algorithm used, but the direction itself is not physically
# relevant, so we get the signs to flip any anti-parallel vectors before
# comparing the two results arrays
signs = np.sign(np.einsum("ij,ij->i", result, ref))
assert_almost_equal(result * signs[:, np.newaxis], ref)

@pytest.fixture()
def universe_pa(self):
Expand Down
Loading