Skip to content

Commit 51af998

Browse files
authored
fix(base): angle_std returned nan for identical angles (#254)
## The problem `angle_std` of a set of identical angles can return `nan`: ```python >>> import numpy as np >>> from spatialmath.base import angle_std >>> angle_std(np.full(100, -2.9)) nan # RuntimeWarning: invalid value encountered in sqrt ``` `R`, the length of the mean of the unit vectors, is at most 1, but rounding can make it `1 + 2e-16`. Then `-2 * log(R)` is a tiny negative number and its square root is `nan`. The answer should be 0. I counted 33 failures in 378 combinations of angle (`-3.1` to `3.1`) and sample size (1 to 100). ## The fix Clamp the value under the square root at zero. `R == 0` still gives `inf`, as before. ## Test `test_angle_std_identical_angles` runs the 378 combinations and fails on the old code and passes now. The existing `test_angle_stats` still passes (24 tests in the file). Found while reviewing a particle filter change in the Robotics Toolbox that uses `angle_std`.
1 parent 0bb35fa commit 51af998

2 files changed

Lines changed: 16 additions & 1 deletion

File tree

‎spatialmath/base/vectors.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -698,7 +698,10 @@ def angle_std(theta: ArrayLike) -> float:
698698
Y = np.sin(theta).mean()
699699
R = np.sqrt(X**2 + Y**2)
700700

701-
return np.sqrt(-2 * np.log(R))
701+
# R is the length of an average of unit vectors so is at most 1, but
702+
# rounding can make it 1 + eps, which would give the square root of a
703+
# (tiny) negative number and nan
704+
return np.sqrt(max(0.0, -2 * np.log(R)))
702705

703706

704707
def angle_mean(theta: ArrayLike) -> float:

‎tests/base/test_vectors.py‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
88
"""
99

10+
import warnings
1011
import numpy as np
1112
import numpy.testing as nt
1213
import unittest
@@ -397,6 +398,17 @@ def test_angle_stats(self):
397398
self.assertAlmostEqual(angle_mean(theta), pi)
398399
self.assertAlmostEqual(angle_std(theta), 0.9717284050981313)
399400

401+
def test_angle_std_identical_angles(self):
402+
# the length of the mean unit vector can come out as 1 + eps, which
403+
# used to give the square root of a negative number
404+
for angle in np.linspace(-3.1, 3.1, 63):
405+
for n in (1, 2, 3, 7, 10, 100):
406+
with warnings.catch_warnings():
407+
warnings.simplefilter("error")
408+
std = angle_std(np.full(n, angle))
409+
self.assertTrue(np.isfinite(std))
410+
self.assertLess(std, 1e-6)
411+
400412
def test_removesmall(self):
401413
v = np.r_[1, 2, 3]
402414
nt.assert_array_almost_equal(removesmall(v), v)

0 commit comments

Comments
 (0)