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.

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

Jump to bottom

DM-9765: Suspicious numerical precision code in Angle #237

Merged
merged 3 commits into from
May 25, 2017
Merged

DM-9765: Suspicious numerical precision code in Angle #237

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1a80cdc
Fix limit handling in Angle.wrapNear
r-owen May 24, 2017
97a02d8
Modify wrapNear limits handling code
r-owen May 24, 2017
da074f1
Modernized testAngle.py
r-owen May 24, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
15 changes: 7 additions & 8 deletions include/lsst/afw/geom/Angle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#if !defined(LSST_AFW_GEOM_ANGLE_H)
#define LSST_AFW_GEOM_ANGLE_H

#include <limits>
#include <iostream>
#include <type_traits>

Expand Down Expand Up @@ -180,17 +179,16 @@ class Angle final {
Angle wrapCtr() const noexcept;

/**
* Wrap this angle to a value `x` such that -&pi; &le; `x - refAng` < &pi;.
* Wrap this angle to a value `x` such that -&pi; &le; `x - refAng` &le; &pi;, approximately.
*
* @param refAng reference angle to match
*
* @returns an angle in the custom normalized interval.
*
* @exceptsafe Shall not throw exceptions.
*
* @warning Exact limits are only guaranteed for radians; limits for other
* units may be slightly squishy due to roundoff errors. There are known
* violations that are demonstrated in testWrap in tests/angle.py.
* @warning The only firm limit is -&pi; &le; `x - refAng` in radians. The upper limit in radians
* and both limits in other units are somewhat squishy, due to roundoff error.
*/
Angle wrapNear(Angle const& refAng) const noexcept;

Expand Down Expand Up @@ -403,13 +401,14 @@ inline Angle Angle::wrapNear(Angle const& refAng) const noexcept {
double const refAngRad = refAng.asRadians();
double wrapped = (*this - refAng).wrapCtr().asRadians() + refAngRad;

// roundoff can cause slightly out-of-range values; fix those
// roundoff can cause slightly out-of-range values; fix those as bast we can;
// note that both conditionals are wanted, since either or both may execute
// (though if/else could be used if the lower limit was squishy for radians)
if (wrapped - refAngRad >= PI) {
wrapped -= TWOPI;
}
// maximum relative roundoff error for subtraction is 2 epsilon
if (wrapped - refAngRad < -PI) {
wrapped -= wrapped * 2.0 * std::numeric_limits<double>::epsilon();
wrapped += TWOPI;
Copy link
Member

Choose a reason for hiding this comment

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

While this is certainly cleaner, I do wonder why the original code didn't use TWOPI here when it did everywhere else.

Copy link
Contributor Author

@r-owen r-owen May 25, 2017 •
edited

Choose a reason for hiding this comment

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

Because the original code tried to keep within bounds at both ends in radians. By adding 2 pi the result may put wrapped - refAngRad >= 2 pi. Still, I'm inclined to live with that. The limits are squishy for other units anyway and we have no use case for precise limits. I will update the docs for wrapNear. Also the test needs some work to exercise cases where the upper limit is violated.

}
return wrapped * radians;
}
Expand Down
236 changes: 111 additions & 125 deletions tests/testAngle.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
>>> import angle; angle.run()
"""
from __future__ import absolute_import, division, print_function
import itertools
import math
import unittest

Expand Down Expand Up @@ -100,29 +101,22 @@ def testComparison(self):
a2 = 2.0 * afwGeom.arcseconds
a1 = 0.5 * afwGeom.arcseconds
a3 = 0.5 * afwGeom.arcseconds
print('a1', a1)
print('a2', a2)
print('a3', a3)
self.assertEqual(a1 == a3, True)
self.assertEqual(a1 != a2, True)
self.assertEqual(a1 <= a2, True)
self.assertEqual(a1 < a2, True)
self.assertEqual(a2 > a1, True)
self.assertEqual(a2 >= a1, True)

self.assertEqual(a1 != a3, False)
self.assertEqual(a1 == a2, False)
self.assertEqual(a1 >= a2, False)
self.assertEqual(a1 > a2, False)
self.assertEqual(a2 < a1, False)
self.assertEqual(a2 <= a1, False)

self.assertEqual(a1 == None, False)
self.assertEqual(None == a1, False)
self.assertEqual(a1 != None, True)
self.assertEqual(None != a1, True)
self.assertEqual(a1 == float(a1), True)
self.assertEqual(float(a1) == a1, True)
self.assertEqual(a1, a3)
self.assertNotEqual(a1, a2)
self.assertLessEqual(a1, a2)
self.assertLess(a1, a2)
self.assertGreater(a2, a1)
self.assertGreaterEqual(a2, a1)

self.assertFalse(a1 != a3)
self.assertFalse(a1 == a2)
self.assertFalse(a1 >= a2)
self.assertFalse(a1 > a2)
self.assertFalse(a2 < a1)
self.assertFalse(a2 <= a1)

self.assertTrue(a1 == float(a1))
self.assertTrue(float(a1) == a1)

def testTrig(self):
self.assertEqual(math.cos(self.d), -1.0)
Expand Down Expand Up @@ -169,109 +163,101 @@ def checkWrappedAngle(self, observed, expected):

def testWrap(self):
eps = np.finfo(float).eps
oneEightyWithSlop = 180 * (1 + eps)
self.assertNotEqual(1 + eps, eps)
for wrap in (-1000, -10, -1, 0, 1, 10, 1000):
for offset in (-2*math.pi, -math.pi, -math.pi*0.5, 0.0,
math.pi*0.5, math.pi*0.75, math.pi, math.pi*2.0):
for epsMult in (-3, -2, -1, 0, 1, 2, 3):
angRad = (offset + (wrap*math.pi)) * (1 + (eps*epsMult))
ang = angRad * afwGeom.radians
angDeg = ang.asDegrees()
sinAng = math.sin(angRad)
cosAng = math.cos(angRad)

posAng = (angRad * afwGeom.radians).wrap()
posAngRad = posAng.asRadians()
posAngDeg = posAng.asDegrees()
posAngArcmin = posAng.asArcminutes()
posAngArcsec = posAng.asArcseconds()
# the code promises 0 <= posAng for all units
self.assertGreaterEqual(posAngRad, 0)
self.assertGreaterEqual(posAngDeg, 0)
self.assertGreaterEqual(posAngArcmin, 0)
self.assertGreaterEqual(posAngArcsec, 0)
# wrap promises posAng < 2*pi only for radians,
# but it seems to work for all units
self.assertLess(posAngRad, 2*math.pi)
self.assertLess(posAngDeg, 360)
self.assertLess(posAngArcmin, 360 * 60)
self.assertLess(posAngArcsec, 360 * 3600)
# prove that posAngDeg and angDeg are the same angle
posErrAng = ((posAngDeg - angDeg) * afwGeom.degrees) \
.wrapCtr()
self.assertAlmostEqual(posErrAng.asDegrees(), 0)
# a sanity check in case wrapCtr gives the wrong answer
self.assertAlmostEqual(math.sin(posAngRad), sinAng)
self.assertAlmostEqual(math.cos(posAngRad), cosAng)

ctrAng = (angRad * afwGeom.radians).wrapCtr()
ctrAngRad = ctrAng.asRadians()
ctrAngDeg = ctrAng.asDegrees()
ctrAngArcmin = ctrAng.asArcminutes()
ctrAngArcsec = ctrAng.asArcseconds()
# wrapCtr promises -pi <= ctrAngRad < pi only for radians,
# but it seems to work for all units
self.assertGreaterEqual(ctrAngRad, -math.pi)
self.assertGreaterEqual(ctrAngDeg, -180)
self.assertGreaterEqual(ctrAngArcmin, -180 * 60)
self.assertGreaterEqual(ctrAngArcsec, -180 * 3600)
self.assertLess(ctrAngRad, math.pi)
self.assertLess(ctrAngDeg, 180)
self.assertLess(ctrAngArcmin, 180 * 60)
self.assertLess(ctrAngArcsec, 180 * 3600)
# prove that ctrAngDeg and ang are the same angle
ctrErrAng = ((ctrAngDeg - angDeg) * afwGeom.degrees) \
.wrapCtr()
self.assertAlmostEqual(ctrErrAng.asDegrees(), 0)
self.assertAlmostEqual(math.sin(ctrAngRad), sinAng)
self.assertAlmostEqual(math.cos(ctrAngRad), cosAng)

for refAngBase in (-math.pi, 0.0, math.pi, math.pi*2.0):
for refEpsMult in (-3, -2, -1, 0, 1, 2, 3):
refAngRad = refAngBase * (1 + (eps * refEpsMult))
refAng = refAngRad * afwGeom.radians
refAngDeg = refAng.asDegrees()
refAngArcmin = refAng.asArcminutes()
refAngArcsec = refAng.asArcseconds()
nearAng = (angRad * afwGeom.radians) \
.wrapNear(refAng)
nearAngRad = nearAng.asRadians()
nearAngDeg = nearAng.asDegrees()
nearAngArcmin = nearAng.asArcminutes()
nearAngArcsec = nearAng.asArcseconds()
# wrapNear promises nearAngRad - refAngRad >= -pi
# for radians but has known failures due to
# roundoff error for other units
self.assertGreaterEqual(nearAngRad - refAngRad,
-math.pi)
self.assertGreaterEqual(nearAngDeg - refAngDeg,
-oneEightyWithSlop)
self.assertGreaterEqual(
nearAngArcmin - refAngArcmin,
-oneEightyWithSlop * 60)
self.assertGreaterEqual(
nearAngArcsec - refAngArcsec,
-oneEightyWithSlop * 3600)
# wrapNear promises nearAngRad - refAngRad < pi
# for radians but has known failures due to
# roundoff error for other units
self.assertLess(nearAngRad - refAngRad,
math.pi)
self.assertLess(nearAngDeg - refAngDeg,
oneEightyWithSlop)
self.assertLess(nearAngArcmin - refAngArcmin,
oneEightyWithSlop * 60)
self.assertLess(nearAngArcsec - refAngArcsec,
oneEightyWithSlop * 3600)
# prove that nearAng and ang are the same angle
nearErrAng = ((nearAngRad - angRad) *
afwGeom.radians).wrapCtr()
self.assertAlmostEqual(nearErrAng.asRadians(), 0)
self.assertAlmostEqual(math.sin(nearAngRad),
sinAng)
self.assertAlmostEqual(math.cos(nearAngRad),
cosAng)
for wrap, offset, epsMult in itertools.product(
(-1000, -10, -1, 0, 1, 10, 1000),
(-2*math.pi, -math.pi, -math.pi*0.5, 0.0, math.pi*0.5, math.pi*0.75, math.pi, math.pi*2.0),
(-3, -2, -1, 0, 1, 2, 3),
):
angRad = (offset + (wrap*math.pi)) * (1 + (eps*epsMult))
ang = angRad * afwGeom.radians
angDeg = ang.asDegrees()
sinAng = math.sin(angRad)
cosAng = math.cos(angRad)

posAng = (angRad * afwGeom.radians).wrap()
posAngRad = posAng.asRadians()
posAngDeg = posAng.asDegrees()
posAngArcmin = posAng.asArcminutes()
posAngArcsec = posAng.asArcseconds()
# the code promises 0 <= posAng for all units
self.assertGreaterEqual(posAngRad, 0)
self.assertGreaterEqual(posAngDeg, 0)
self.assertGreaterEqual(posAngArcmin, 0)
self.assertGreaterEqual(posAngArcsec, 0)
# wrap promises posAng < 2*pi only for radians,
# but it seems to work for all units
self.assertLess(posAngRad, 2*math.pi)
self.assertLess(posAngDeg, 360)
self.assertLess(posAngArcmin, 360 * 60)
self.assertLess(posAngArcsec, 360 * 3600)
# prove that posAngDeg and angDeg are the same angle
posErrAng = ((posAngDeg - angDeg) * afwGeom.degrees) \
.wrapCtr()
self.assertAlmostEqual(posErrAng.asDegrees(), 0)
# a sanity check in case wrapCtr gives the wrong answer
self.assertAlmostEqual(math.sin(posAngRad), sinAng)
self.assertAlmostEqual(math.cos(posAngRad), cosAng)

ctrAng = (angRad * afwGeom.radians).wrapCtr()
ctrAngRad = ctrAng.asRadians()
ctrAngDeg = ctrAng.asDegrees()
ctrAngArcmin = ctrAng.asArcminutes()
ctrAngArcsec = ctrAng.asArcseconds()
# wrapCtr promises -pi <= ctrAngRad < pi only for radians,
# but it seems to work for all units
self.assertGreaterEqual(ctrAngRad, -math.pi)
self.assertGreaterEqual(ctrAngDeg, -180)
self.assertGreaterEqual(ctrAngArcmin, -180 * 60)
self.assertGreaterEqual(ctrAngArcsec, -180 * 3600)
self.assertLess(ctrAngRad, math.pi)
self.assertLess(ctrAngDeg, 180)
self.assertLess(ctrAngArcmin, 180 * 60)
self.assertLess(ctrAngArcsec, 180 * 3600)
# prove that ctrAngDeg and ang are the same angle
ctrErrAng = ((ctrAngDeg - angDeg) * afwGeom.degrees) \
.wrapCtr()
self.assertAlmostEqual(ctrErrAng.asDegrees(), 0)
self.assertAlmostEqual(math.sin(ctrAngRad), sinAng)
self.assertAlmostEqual(math.cos(ctrAngRad), cosAng)

for refAngBase, refAngWrap, refEpsMult in itertools.product(
(-math.pi, 0.0, math.pi, math.pi*2.0),
range(-10, 11, 2),
(-3, -2, -1, 0, 1, 2, 3),
):
refAngRad = refAngBase * (1 + (eps * refEpsMult)) + refAngWrap * math.pi * 2.0
refAng = refAngRad * afwGeom.radians
refAngDeg = refAng.asDegrees()
refAngArcmin = refAng.asArcminutes()
refAngArcsec = refAng.asArcseconds()
nearAng = (angRad * afwGeom.radians).wrapNear(refAng)
nearAngRad = nearAng.asRadians()
nearAngDeg = nearAng.asDegrees()
nearAngArcmin = nearAng.asArcminutes()
nearAngArcsec = nearAng.asArcseconds()

fudgeFactor = 5 # 1 and 2 are too small for one case, 3 works; 5 has margin
relEps = max(1, nearAngRad / math.pi, refAngRad / math.pi) * eps * fudgeFactor
piWithSlop = math.pi * (1 + relEps)
oneEightyWithSlop = 180 * (1 + relEps)

# wrapNear promises nearAngRad - refAngRad >= -pi
# for radians but has known failures due to
# roundoff error for other units and for < pi
self.assertGreaterEqual(nearAngRad - refAngRad, -math.pi)
self.assertGreaterEqual(nearAngDeg - refAngDeg, -oneEightyWithSlop)
self.assertGreaterEqual(nearAngArcmin - refAngArcmin, -oneEightyWithSlop * 60)
self.assertGreaterEqual(nearAngArcsec - refAngArcsec, -oneEightyWithSlop * 3600)
self.assertLessEqual(nearAngRad - refAngRad, piWithSlop)
self.assertLessEqual(nearAngDeg - refAngDeg, oneEightyWithSlop)
self.assertLessEqual(nearAngArcmin - refAngArcmin, oneEightyWithSlop * 60)
self.assertLessEqual(nearAngArcsec - refAngArcsec, oneEightyWithSlop * 3600)
# prove that nearAng and ang are the same angle
nearErrAng = ((nearAngRad - angRad) * afwGeom.radians).wrapCtr()
self.assertAlmostEqual(nearErrAng.asRadians(), 0)
self.assertAlmostEqual(math.sin(nearAngRad), sinAng)
self.assertAlmostEqual(math.cos(nearAngRad), cosAng)


class MemoryTester(lsst.utils.tests.MemoryTestCase):
Expand Down