Float comparison issues with time and quantity

The following snippet returns True with astropy 2.0.2, but False on master and 2.0.x (2.0.3.dev19469)

>>> import astropy
>>> print(astropy.__version__)
>>> from astropy.time import Time
>>> import astropy.units as u
>>> start = Time('2016-02-06 03:00:00')
>>> duration = 2*u.hour
>>> end_time = start + duration
>>> end_time - duration == start
2.0.2
True

>>> import astropy
>>> print(astropy.__version__)
>>> from astropy.time import Time
>>> import astropy.units as u
>>> start = Time('2016-02-06 03:00:00')
>>> duration = 2*u.hour
>>> end_time = start + duration
>>> print(end_time - duration == start)
>>> print(end_time - duration)
>>> print(start)
3.0.dev20874
False
2016-02-06 03:00:00.000
2016-02-06 03:00:00.000

This is weird. Same problem on master, but

(end_time - duration - start).to(u.ns)
# <Quantity 0.0 ns>

It clearly is some kind of rounding problem, though:

In [11]: start.jd1, start.jd2
Out[11]: (2457425.0, -0.375)

In [12]: (end_time - duration).jd1, (end_time - duration).jd2
Out[12]: (2457425.0, -0.37500000000000006)

Actually, it is a bit of a stranger rounding problem: since hours are not a constant fraction of a day, the time is converted to TAI before the duration is subtracted, and then converted back to UTC. The problem doesn't happen if we start with TAI:

In [15]: start = Time('2016-02-06 03:00:00', scale='tai')

In [16]: end_time = start + duration

In [17]: end_time - duration == start
Out[17]: True

I think the underlying cause of the change is #6653 - I think the problem it fixed is greater than this one, though.

p.s. integer for Time is days; 2 hours = 1/12 of a day, so not exactly representable in binary.

A bit further investigation: it is a round-trip error in converting from utc to tai:

start = Time('2016-02-06 03:00:00')
start == start.tai.utc
# False
start.jd2
# -0.375
start.tai.utc.jd2
# -0.37500000000000006

Note that it is just unlucky rounding. The same happens in 2.0.2 if I ensure input is rounded as it normally is:

start = Time(2457424.625, format='jd')
start.jd1, start.jd2
# (2457425.0, -0.375)
start == start.tai.utc
# False
start = Time('2016-02-06 03:00:00')
start.jd1, start.jd2
# (2457424.5, 0.125)
start == start.tai.utc
# True

Checking even further, it boils down to erfa not quite round-tripping:

import astropy._erfa as erfa
erfa.taiutc(*erfa.utctai(2457424.5, 0.125))
# (array(2457424.5), array(0.125))
erfa.taiutc(*erfa.utctai(2457425.0, -.375))
# (array(2457425.0), array(-0.37500000000000006))

I can make the particular problem go away with a small change in jd2cal. Arguably, it is now an upstream problem.

@mhvk

@mhvk - Thanks for investigating. Take away for downstream will be to adjust those tests as all of it is upstream from astroplan.

Yes, I think any tests like this on time should test closeness (<1ns is reasonable) rather than equality.

Yes, I think any tests like this on time should test closeness (<1ns is reasonable) rather than equality.

Actually I'm not sure I agree there. You never know when someone really wanted to check that it's closer than 1 ns...? Same reason why we don't do equality in coordinates et al - "equals" really is ambiguous here if you don't mean literally bit-wise equality...?

@eteq

@eteq - the comment was about tests on round-tripping. Since this depends on float rounding working the way we want, I don't think we can guarantee equality, so all I am suggesting is that we do not use __eq__ in tests, not that we adjust the meaning of __eq__!

Maybe it's because of #6653?

@mhvk

Sorry for the noise, @mhvk said exactly that - should have "Ctrl+F"-ed.

Looking at older issues, I think the summary here is whether it is worth trying to improve round-tripping for ERFA/SOFA, and push that upstream.

Sent the following to SOFA:

Hi Patrick,

I understand you have been in touch last Summer with Tim Jenness (CC'd)
about a small round-trip problem in astropy, which we traced back to
SOFA (see
https://github.com/astropy/astropy/issues/6970#issuecomment-351461606).

I found a fix for it already at the time, and vaguely recall trying to
pass it on, but apparently it never made it, so let me try here again.

The problem is that UTC->TAI->UTC does not roundtrip, even though there
really is nothing to prevent it from doing so if one is careful with how
floats are dealt with.  In this case, the patch is very simple:

diff jd2cal.c~ jd2cal.c
85c85
<    d2 -= 0.5;
---
>    d1 -= 0.5;

Here, the trick is to subtract 0.5 from the larger of the two doubles.
Since this is always a large positive number, the subtraction of 0.5 is
guaranteed to be exact, while for the smaller of the two doublds, this
is not the case.  The same change should be made in jdcalf.c.:

diff jdcalf.c~ jdcalf.c
91c91
<    d2 -= 0.5;
---
>    d1 -= 0.5;

For a regression test, you might want to use the attached additional
test (note the negative input value of d2; by subtracting 0.5 one looses
one bit of precision).  Note that this is for the specific issue we
found, so it covers jd2cal (indirectly) but not jdcalf.  It might be
better to have a specific test for jd2cal itself.

With all best wishes,

Marten

p.s. It would be so good if SOFA got an issue tracker, or even better
     allowed pull requests, or at least have a git version so one can
     send e-mail patches.

--
Prof. M. H. van Kerkwijk
Dept. of Astronomy & Astroph., 50 St George St., Toronto, ON, M5S 3H4, Canada
McLennan Labs, room 1203B, tel: +1(416)9467288, fax: +1(416)9467287
[email protected], http://www.astro.utoronto.ca/~mhvk


diff t_sofa_c.c~ t_sofa_c.c
9630a9631,9661
> static void t_utctaiutc(int *status)
> /*
> **  - - - - - - - - -
> **   t _ u t c t a i
> **  - - - - - - - - -
> **
> **  Test iauUtctai and iauTaiutc functions, verifying roundtrip.
> **
> **  Returned:
> **     status    int         FALSE = success, TRUE = fail
> **
> **  Called:  iauUtctai, iauTaiutc, vvd, viv
> **
> **  This revision:  2020 April 8
> */
> {
>    double u1, u2, t1, t2, ru1, ru2;
>    int j;
>
>    u1 = 2457425.0;
>    u2 = -0.375;
>
>    j = iauUtctai(u1, u2, &t1, &t2);
>    viv(j, 0, "iauUtctai", "j", status);
>    j = iauTaiutc(t1, t2, &ru1, &ru2);
>    viv(j, 0, "iauTaiutc", "j", status);
>    vvd(ru1, u1, 0., "iauUtctai,iauTaiutc", "ru1", status);
>    vvd(ru2, u2, 0., "iauUtctai,iauTaiutc", "ru2", status);
>
> }
>
10140a10172
>    t_utctaiutc(&status);

After quite a bit of back and forth with Pat, he convinced me that perfect roundtripping was simply impossible, since adding and subtracting the UTC-TAI difference cannot always be done without precision loss. But with fixes to jd2cal (the above plus some extra), at least the maximum precision loss has been reduced by a factor 4.

Given this, I think the remaining to-do issue is to add documentation that for times in UTC, round-tripping errors can cause even adding 0 seconds to give a non-identical time.

Navigation Menu

Repository navigation

Float comparison issues with time and quantity #6970

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Float comparison issues with time and quantity #6970

Description

Activity

mhvk commented on Dec 12, 2017

mhvk commented on Dec 12, 2017

mhvk commented on Dec 12, 2017

mhvk commented on Dec 12, 2017

mhvk commented on Dec 13, 2017

bsipocz commented on Dec 13, 2017

mhvk commented on Dec 13, 2017

eteq commented on Mar 15, 2018

mhvk commented on Mar 15, 2018

astrojuanlu commented on May 13, 2018

astrojuanlu commented on May 13, 2018

mhvk commented on Mar 25, 2019

mhvk commented on Apr 8, 2020

mhvk commented on Apr 10, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions