Floating-point-related failures on arm64/aarch64 #9661
Description
Description
On arm64/aarch64 (essentially alternative names for the same ABI, one of two implemented on the 64-bit ARM CPUs in modern phones and Raspberry PIs) double precision calculations do not have access to the hidden extra 16 bits that they do on Intel machines, and long double is actually quadruple precision, alas implemented in software. As a result a number of floating-point differences arise, chiefly when using Time to push the limits of available precision. (In fact np.longdouble here has more precision than Time.) PR #9407 will add arm64 tests to CI but unfortunately amr64 support on Travis and on PyPI is limited (no compiled wheels for anything, no conda) so that can't be implemented yet. But some arm64 machines are affordable, so I ran our test suite and uncovered a number of problems, mostly in the test suite.
A few of the errors below are our usual test failures when run with remote data, a few are tests that assume long doubles are less precise than Time, a few result from different float error handling, but I've included everything because it wasn't clear to me which was which or whether there was another category.
Expected behavior
Test suite passes.
Actual behavior
====================================================== FAILURES =======================================================
______________________________________ test_database_specify[NGC 3642-db_dict0] _______________________________________
name = 'NGC 3642'
db_dict = {'all': '# ngc3642 #Q22523722\n#=S=Simbad (via url): 1\n%@ 503952\n%I.0 NGC 3642\n%C.0 LIN\n%C.N0 15.15
.01.00\n%...eR (local): 1\n%J 170.56 +59.08 = 11:22.2 +59:05\n%I.0 {NGC} 3642\n\n\n\n#====Done (2013-Feb-12,16:3
7:42z)===='}
@pytest.mark.remote_data
@pytest.mark.parametrize(("name", "db_dict"), [('NGC 3642', _cached_ngc3642),
('castor', _cached_castor)])
def test_database_specify(name, db_dict):
# First check that at least some sesame mirror is up
for url in sesame_url.get():
if urllib.request.urlopen(url).getcode() == 200:
break
else:
pytest.skip("All SESAME mirrors appear to be down, skipping "
"test_name_resolve.py:test_database_specify()...")
for db in db_dict.keys():
with sesame_database.set(db):
> icrs = SkyCoord.from_name(name)
astropy/coordinates/tests/test_name_resolve.py:167:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
astropy/coordinates/sky_coordinate.py:1697: in from_name
icrs_coord = get_icrs_coordinates(name, parse)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
name = 'NGC 3642', parse = False
def get_icrs_coordinates(name, parse=False):
"""
Retrieve an ICRS object by using an online name resolving service to
retrieve coordinates for the specified name. By default, this will
search all available databases until a match is found. If you would like
to specify the database, use the science state
``astropy.coordinates.name_resolve.sesame_database``. You can also
specify a list of servers to use for querying Sesame using the science
state ``astropy.coordinates.name_resolve.sesame_url``. This will try
each one in order until a valid response is returned. By default, this
list includes the main Sesame host and a mirror at vizier. The
configuration item `astropy.utils.data.Conf.remote_timeout` controls the
number of seconds to wait for a response from the server before giving
up.
Parameters
----------
name : str
The name of the object to get coordinates for, e.g. ``'M42'``.
parse: bool
Whether to attempt extracting the coordinates from the name by
parsing with a regex. For objects catalog names that have
J-coordinates embedded in their names eg:
'CRTS SSS100805 J194428-420209', this may be much faster than a
sesame query for the same object name. The coordinates extracted
in this way may differ from the database coordinates by a few
deci-arcseconds, so only use this option if you do not need
sub-arcsecond accuracy for coordinates.
Returns
-------
coord : `astropy.coordinates.ICRS` object
The object's coordinates in the ICRS frame.
"""
# if requested, first try extract coordinates embedded in the object name.
# Do this first since it may be much faster than doing the sesame query
if parse:
from . import jparser
if jparser.search(name):
return jparser.to_skycoord(name)
else:
# if the parser failed, fall back to sesame query.
pass
# maybe emit a warning instead of silently falling back to sesame?
database = sesame_database.get()
# The web API just takes the first letter of the database name
db = database.upper()[0]
# Make sure we don't have duplicates in the url list
urls = []
domains = []
for url in sesame_url.get():
domain = urllib.parse.urlparse(url).netloc
# Check for duplicates
if domain not in domains:
domains.append(domain)
# Add the query to the end of the url, add to url list
fmt_url = os.path.join(url, "{db}?{name}")
fmt_url = fmt_url.format(name=urllib.parse.quote(name), db=db)
urls.append(fmt_url)
exceptions = []
for url in urls:
try:
# Retrieve ascii name resolve data from CDS
resp = urllib.request.urlopen(url, timeout=data.conf.remote_timeout)
resp_data = resp.read()
break
except urllib.error.URLError as e:
exceptions.append(e)
continue
except socket.timeout as e:
# There are some cases where urllib2 does not catch socket.timeout
# especially while receiving response data on an already previously
# working request
e.reason = "Request took longer than the allowed {:.1f} " \
"seconds".format(data.conf.remote_timeout)
exceptions.append(e)
continue
# All Sesame URL's failed...
else:
messages = [f"{url}: {e.reason}"
for url, e in zip(urls, exceptions)]
raise NameResolveError("All Sesame queries failed. Unable to "
"retrieve coordinates. See errors per URL "
"below: \n {}".format("\n".join(messages)))
ra, dec = _parse_response(resp_data)
if ra is None and dec is None:
if db == "A":
err = f"Unable to find coordinates for name '{name}'"
else:
err = "Unable to find coordinates for name '{}' in database {}"\
.format(name, database)
> raise NameResolveError(err)
E astropy.coordinates.name_resolve.NameResolveError: Unable to find coordinates for name 'NGC 3642' in databa
se vizier
astropy/coordinates/name_resolve.py:191: NameResolveError
___________________________ [doctest] astropy.stats.info_theory.bayesian_info_criterion_lsq ___________________________
171 >>> g_init = models.Gaussian1D(amplitude=1., mean=0, stddev=1.)
172 >>> fit_g = fitting.LevMarLSQFitter()
173 >>> g = fit_g(g_init, x, y)
174 >>> # Compute the mean squared errors
175 >>> ssr_t = np.sum((t(x) - y)*(t(x) - y))
176 >>> ssr_g = np.sum((g(x) - y)*(g(x) - y))
177 >>> # Compute the bics
178 >>> bic_t = bayesian_info_criterion_lsq(ssr_t, 4, x.shape[0])
179 >>> bic_g = bayesian_info_criterion_lsq(ssr_g, 3, x.shape[0])
180 >>> bic_t - bic_g # doctest: +FLOAT_CMP
Expected:
30.644474706065466
Got:
31.11155402550662
/tmp/astropy-test-0eiyyonh/lib/python3.7/site-packages/astropy/stats/info_theory.py:180: DocTestFailure
____________________________________ TestNumericalSubFormat.test_explicit_example _____________________________________
self = <astropy.time.tests.test_basic.TestNumericalSubFormat object at 0xffff79c23748>
def test_explicit_example(self):
t = Time('54321.000000000001', format='mjd')
assert t == Time(54321, 1e-12, format='mjd')
assert t.mjd == 54321. # Lost precision!
assert t.value == 54321. # Lost precision!
assert t.to_value('mjd') == 54321. # Lost precision!
assert t.to_value('mjd', subfmt='str') == '54321.000000000001'
assert t.to_value('mjd', 'bytes') == b'54321.000000000001'
expected_long = np.longdouble(54321.) + np.longdouble(1e-12)
> assert t.to_value('mjd', subfmt='long') == expected_long
E AssertionError: assert 54321.000000000000999977878279878496 == 54321.00000000000099999999999999998
E + where 54321.000000000000999977878279878496 = <bound method Time.to_value of <Time object: scale='utc' forma
t='mjd' value=54321.0>>('mjd', subfmt='long')
E + where <bound method Time.to_value of <Time object: scale='utc' format='mjd' value=54321.0>> = <Time objec
t: scale='utc' format='mjd' value=54321.0>.to_value
astropy/time/tests/test_basic.py:870: AssertionError
___________________________________ TestNumericalSubFormat.test_explicit_longdouble ___________________________________
self = <astropy.time.tests.test_basic.TestNumericalSubFormat object at 0xffff79e92da0>
@pytest.mark.skipif(np.finfo(np.longdouble).eps >= np.finfo(float).eps,
reason="long double is the same as float")
def test_explicit_longdouble(self):
i = 54321
f = 2.**(-np.finfo(np.longdouble).nmant) * 65536
mjd_long = np.longdouble(i) + np.longdouble(f)
assert mjd_long != i, "longdouble failure!"
t = Time(mjd_long, format='mjd')
expected = Time(i, f, format='mjd')
assert t == expected
t_float = Time(i+f, format='mjd')
assert t_float == Time(i, format='mjd')
> assert t_float != t
E AssertionError: assert <Time object: scale='utc' format='mjd' value=54321.0> != <Time object: scale='utc' forma
t='mjd' value=54321.0>
astropy/time/tests/test_basic.py:896: AssertionError
_____________________________ TestNumericalSubFormat.test_longdouble_for_other_types[mjd] _____________________________
self = <astropy.time.tests.test_basic.TestNumericalSubFormat object at 0xffff7d0bbef0>, fmt = 'mjd'
@pytest.mark.skipif(np.finfo(np.longdouble).eps >= np.finfo(float).eps,
reason="long double is the same as float")
@pytest.mark.parametrize("fmt", ["mjd", "unix", "cxcsec"])
def test_longdouble_for_other_types(self, fmt):
t_fmt = getattr(Time(58000, format="mjd"), fmt) # Get regular float
t_fmt_long = np.longdouble(t_fmt)
t_fmt_long2 = t_fmt_long * (np.finfo(np.longdouble).eps * 2 + 1)
assert t_fmt_long != t_fmt_long2, "longdouble weird!"
tm = Time(t_fmt_long, format=fmt)
tm2 = Time(t_fmt_long2, format=fmt)
> assert tm != tm2
E AssertionError: assert <Time object: scale='utc' format='mjd' value=58000.0> != <Time object: scale='utc' forma
t='mjd' value=58000.0>
astropy/time/tests/test_basic.py:912: AssertionError
____________________________ TestNumericalSubFormat.test_longdouble_for_other_types[unix] _____________________________
self = <astropy.time.tests.test_basic.TestNumericalSubFormat object at 0xffff79c306d8>, fmt = 'unix'
@pytest.mark.skipif(np.finfo(np.longdouble).eps >= np.finfo(float).eps,
reason="long double is the same as float")
@pytest.mark.parametrize("fmt", ["mjd", "unix", "cxcsec"])
def test_longdouble_for_other_types(self, fmt):
t_fmt = getattr(Time(58000, format="mjd"), fmt) # Get regular float
t_fmt_long = np.longdouble(t_fmt)
t_fmt_long2 = t_fmt_long * (np.finfo(np.longdouble).eps * 2 + 1)
assert t_fmt_long != t_fmt_long2, "longdouble weird!"
tm = Time(t_fmt_long, format=fmt)
tm2 = Time(t_fmt_long2, format=fmt)
> assert tm != tm2
E AssertionError: assert <Time object: scale='utc' format='unix' value=1504483200.0> != <Time object: scale='utc'
format='unix' value=1504483200.0>
astropy/time/tests/test_basic.py:912: AssertionError
___________________________ TestNumericalSubFormat.test_longdouble_for_other_types[cxcsec] ____________________________
self = <astropy.time.tests.test_basic.TestNumericalSubFormat object at 0xffff79c235f8>, fmt = 'cxcsec'
@pytest.mark.skipif(np.finfo(np.longdouble).eps >= np.finfo(float).eps,
reason="long double is the same as float")
@pytest.mark.parametrize("fmt", ["mjd", "unix", "cxcsec"])
def test_longdouble_for_other_types(self, fmt):
t_fmt = getattr(Time(58000, format="mjd"), fmt) # Get regular float
t_fmt_long = np.longdouble(t_fmt)
t_fmt_long2 = t_fmt_long * (np.finfo(np.longdouble).eps * 2 + 1)
assert t_fmt_long != t_fmt_long2, "longdouble weird!"
tm = Time(t_fmt_long, format=fmt)
tm2 = Time(t_fmt_long2, format=fmt)
> assert tm != tm2
E AssertionError: assert <Time object: scale='tt' format='cxcsec' value=620870469.184> != <Time object: scale='tt
' format='cxcsec' value=620870469.184>
astropy/time/tests/test_basic.py:912: AssertionError
_______________________________________ test_mjd_longdouble_preserves_precision _______________________________________
custom_format_name = 'custom_format_name'
def test_mjd_longdouble_preserves_precision(custom_format_name):
class CustomMJD(TimeFormat):
name = custom_format_name
def _check_val_type(self, val, val2):
val = np.longdouble(val)
if val2 is not None:
raise ValueError("Only one value permitted")
return val, 0
def set_jds(self, val, val2):
mjd1 = np.float64(np.floor(val))
mjd2 = np.float64(val - mjd1)
self.jd1, self.jd2 = day_frac(mjd1 + DJM0, mjd2)
@property
def value(self):
mjd1, mjd2 = day_frac(self.jd1 - DJM0, self.jd2)
return np.longdouble(mjd1) + np.longdouble(mjd2)
m = 58000.0
t = Time(m, format=custom_format_name)
t2 = Time(m + 2 * m * np.finfo(np.longdouble).eps, format=custom_format_name)
> assert t != t2
E AssertionError: assert <Time object: scale='utc' format='custom_format_name' value=58000.0> != <Time object: sc
ale='utc' format='custom_format_name' value=58000.0>
astropy/time/tests/test_custom_formats.py:159: AssertionError
___________________________________________ TestComparisonUfuncs.test_sign ____________________________________________
self = <astropy.units.tests.test_quantity_ufuncs.TestComparisonUfuncs object at 0xffff78f8d828>
def test_sign(self):
q = [1., np.inf, -np.inf, np.nan, -1., 0.] * u.m
# Ignore "invalid value encountered in sign" warning on Windows.
if sys.platform.startswith('win'):
ctx = np.errstate(invalid='ignore')
else:
ctx = nullcontext()
with ctx:
> out = np.sign(q)
astropy/units/tests/test_quantity_ufuncs.py:723:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <Quantity [ 1., inf, -inf, nan, -1., 0.] m>, function = <ufunc 'sign'>, method = '__call__'
inputs = (<Quantity [ 1., inf, -inf, nan, -1., 0.] m>,), kwargs = {}, converters = [None], unit = None
out = None, arrays = [array([ 1., inf, -inf, nan, -1., 0.])]
input_ = array([ 1., inf, -inf, nan, -1., 0.]), converter = None
def __array_ufunc__(self, function, method, *inputs, **kwargs):
"""Wrap numpy ufuncs, taking care of units.
Parameters
----------
function : callable
ufunc to wrap.
method : str
Ufunc method: ``__call__``, ``at``, ``reduce``, etc.
inputs : tuple
Input arrays.
kwargs : keyword arguments
As passed on, with ``out`` containing possible quantity output.
Returns
-------
result : `~astropy.units.Quantity`
Results of the ufunc, with the unit set properly.
"""
# Determine required conversion functions -- to bring the unit of the
# input to that expected (e.g., radian for np.sin), or to get
# consistent units between two inputs (e.g., in np.add) --
# and the unit of the result (or tuple of units for nout > 1).
converters, unit = converters_and_unit(function, method, *inputs)
out = kwargs.get('out', None)
# Avoid loop back by turning any Quantity output into array views.
if out is not None:
# If pre-allocated output is used, check it is suitable.
# This also returns array view, to ensure we don't loop back.
if function.nout == 1:
out = out[0]
out_array = check_output(out, unit, inputs, function=function)
# Ensure output argument remains a tuple.
kwargs['out'] = (out_array,) if function.nout == 1 else out_array
# Same for inputs, but here also convert if necessary.
arrays = []
for input_, converter in zip(inputs, converters):
input_ = getattr(input_, 'value', input_)
arrays.append(converter(input_) if converter else input_)
# Call our superclass's __array_ufunc__
> result = super().__array_ufunc__(function, method, *arrays, **kwargs)
E RuntimeWarning: invalid value encountered in sign
astropy/units/quantity.py:481: RuntimeWarning
_________________________________________ TestInplaceUfuncs.test_sign_inplace _________________________________________
self = <astropy.units.tests.test_quantity_ufuncs.TestInplaceUfuncs object at 0xffff78fa0898>
def test_sign_inplace(self):
q = [1., np.inf, -np.inf, np.nan, -1., 0.] * u.m
check = np.empty(q.shape, q.dtype)
# Ignore "invalid value encountered in sign" warning on Windows.
if sys.platform.startswith('win'):
ctx = np.errstate(invalid='ignore')
else:
ctx = nullcontext()
with ctx:
> np.sign(q.value, out=check)
E RuntimeWarning: invalid value encountered in sign
astropy/units/tests/test_quantity_ufuncs.py:931: RuntimeWarning
_______________________________________________ test_footprint_contains _______________________________________________
def test_footprint_contains():
"""
Test WCS.footprint_contains(skycoord)
"""
header = """
WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 1045.0 / Pixel coordinate of reference point
CRPIX2 = 1001.0 / Pixel coordinate of reference point
PC1_1 = -0.00556448550786 / Coordinate transformation matrix element
PC1_2 = -0.001042120133257 / Coordinate transformation matrix element
PC2_1 = 0.001181477028705 / Coordinate transformation matrix element
PC2_2 = -0.005590809742987 / Coordinate transformation matrix element
CDELT1 = 1.0 / [deg] Coordinate increment at reference point
CDELT2 = 1.0 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'RA---TAN' / TAN (gnomonic) projection + SIP distortions
CTYPE2 = 'DEC--TAN' / TAN (gnomonic) projection + SIP distortions
CRVAL1 = 250.34971683647 / [deg] Coordinate value at reference point
CRVAL2 = 2.2808772582495 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 2.2808772582495 / [deg] Native latitude of celestial pole
RADESYS = 'ICRS' / Equatorial coordinate system
MJD-OBS = 58612.339199259 / [d] MJD of observation matching DATE-OBS
DATE-OBS= '2019-05-09T08:08:26.816Z' / ISO-8601 observation date matching MJD-OB
NAXIS = 2 / NAXIS
NAXIS1 = 2136 / length of first array dimension
NAXIS2 = 2078 / length of second array dimension
""" # noqa
header = fits.Header.fromstring(header.strip(), '\n')
test_wcs = wcs.WCS(header)
hasCoord = test_wcs.footprint_contains(SkyCoord(254, 2, unit='deg'))
assert hasCoord
hasCoord = test_wcs.footprint_contains(SkyCoord(240, 2, unit='deg'))
assert not hasCoord
# Ignore "invalid value encountered in less" warning on Windows.
if sys.platform.startswith('win'):
ctx = np.errstate(invalid='ignore')
else:
ctx = nullcontext()
with ctx:
> hasCoord = test_wcs.footprint_contains(SkyCoord(24, 2, unit='deg'))
astropy/wcs/tests/test_wcs.py:1227:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
astropy/wcs/wcs.py:3111: in footprint_contains
return coord.contained_by(self, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <SkyCoord (ICRS): (ra, dec) in deg
(24., 2.)>
wcs = WCS Keywords
Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 250.34971683647 2.2808772582495
CRPIX ...0786 -0.001042120133257
PC2_1 PC2_2 : 0.001181477028705 -0.005590809742987
CDELT : 1.0 1.0
NAXIS : 2136 2078
image = None, kwargs = {}, ymax = 2078, xmax = 2136
warnings = <module 'warnings' from '/home/peridot/.virtualenvs/astropy/lib/python3.7/warnings.py'>, x = array(nan)
y = array(nan)
def contained_by(self, wcs, image=None, **kwargs):
"""
Determines if the SkyCoord is contained in the given wcs footprint.
Parameters
----------
wcs : `~astropy.wcs.WCS`
The coordinate to check if it is within the wcs coordinate.
image : array
Optional. The image associated with the wcs object that the cooordinate
is being checked against. If not given the naxis keywords will be used
to determine if the coordinate falls within the wcs footprint.
**kwargs :
Additional arguments to pass to `~astropy.coordinates.SkyCoord.to_pixel`
Returns
-------
response : bool
True means the WCS footprint contains the coordinate, False means it does not.
"""
if image is not None:
ymax, xmax = image.shape
else:
xmax, ymax = wcs._naxis
import warnings
with warnings.catch_warnings():
# Suppress warnings since they just mean we didn't find the coordinate
warnings.simplefilter("ignore")
try:
x, y = self.to_pixel(wcs, **kwargs)
except Exception:
return False
> return (x < xmax) & (x > 0) & (y < ymax) & (y > 0)
E RuntimeWarning: invalid value encountered in less
astropy/coordinates/sky_coordinate.py:1439: RuntimeWarning
======================= 11 failed, 14523 passed, 108 skipped, 52 xfailed in 3187.82s (0:53:07) ========================Steps to Reproduce
$ python setup.py test --remote-data=any
System Details
Python 3.7.3 (default, Apr 3 2019, 05:39:12)
[GCC 8.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import platform; print(platform.platform())
Linux-5.3.11-rockchip64-aarch64-with-debian-10.2
>>> import sys; print("Python", sys.version)
Python 3.7.3 (default, Apr 3 2019, 05:39:12)
[GCC 8.3.0]
>>> import numpy; print("Numpy", numpy.__version__)
Numpy 1.17.4
>>> import scipy; print("Scipy", scipy.__version__)
Scipy 1.3.3
>>> import astropy; print("astropy", astropy.__version__)
astropy 4.1.dev27006