Time equality comparisons are possibly too precise

In astro_metadata_translator I've been adding some tests to make sure that the end date of an observation comes after the start date of an observation. Sometimes DATE-END header is missing so the end date is derived by adding the exposure time to the start date. This is fine until biases are included which have 0.0 second exposures, whereupon the end date is now in the past.

 $ python
Python 3.7.3 (default, Mar 27 2019, 16:54:48) 
[Clang 4.0.1 (tags/RELEASE_401/final)] :: Anaconda, Inc. on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from astropy.time import Time
>>> import astropy.units as u
>>> t1 = Time("2018-07-24T10:41:56.807015240")
>>> t2 = t1 + 0.0*u.s
>>> from astropy.time import TimeDelta
>>> t3 = t1 + TimeDelta(0.0*u.s)
>>> t3 == t1
False
>>> t2 < t1
True
>>> t1.jd2
-0.054203622508796234
>>> t2.jd2
-0.054203622508796345
>>> t3.jd2
-0.054203622508796345

i.e. if you add (at least on my Mac) 0 seconds to a Time the resulting time object is, in this case, 10 picoseconds behind the original time. This breaks my tests because the end date is not allowed to be behind the start date. I'm fixing the tests by converting the Time to a ISO string and then converting back to a Time object before doing the >= test but I'm wondering whether this behavior is deliberate or if no-one has ever tried it.

Should __eq__ be comparing floating point numbers to 0.0 without some form of precision check (math.isclose?).
If a precision is specified for the Time should that be used in the equality?
If a time delta of 0 seconds is added to a Time should the new Time go backwards in time?
Is Time meant to support picoseconds precision? If it isn't then the equality test should take that into account. If it is then the above suggests it's not.

The precision Time supports is the precision of a double in days, so

(np.finfo(np.float64).eps * u.day).to(u.ps)
# <Quantity 19.18465387 ps>

This is what most erfa functions work with, so not easy to improve.

That said, adding 0 should not result in a time change! The problem here may in fact be that to add anything to a UTC time, first the time is converted to TAI, then the addition done, and finally the result converted back to UTC. And, indeed, it is in the conversion that there is a problem:

t1 = Time("2018-07-24T10:41:56.807015240")
t1.tai.utc < t1
# True

Which reminds me of something we encountered earlier... Yes, this is a duplicate of #6970, where the conclusion was that we could solve this issue with a small change in SOFA/ERFA. Indeed, the fix is trivial - see https://github.com/astropy/astropy/compare/master...mhvk:jd2cal-precision-test?expand=1

If anyone has the time to push this upstream... (we're behind, but it is not in the 2018-01-30 release)

p.s. Unlike #6970, I do count this as a bug: t + 0 * u.s == t should always be True!

p.p.s. That said, it would be good to have a version of isclose or so - I end up writing too many tests of the sort of assert abs(t1 - t2) < 1 * u.ns...

we're behind, but it is not in the 2018-01-30 release

should that be opened as an issue here or upstream, or both?

I think it should go direct to SOFA - we're trying not to deviate from it in ERFA.

Oh, what I meant is that presumably we should update erfa in astropy (but first it also needs to be updated to be the latest sofa release)

Yes, we should update erfa first - see liberfa/erfa#50 - though you're doing enough release work already!! (We may still want to put this in the astropy team and perhaps ask for volunteers...)

@scottransom

For the record, I emailed SOFA a couple of days ago to ask them about the proposed tweak. No reply yet but thought I should note it here.

cc/ @scottransom

Following up, no response from SOFA so maybe I'm using a bad email address.

I did get a response from Pat Wallace but unfortunately it came during my summer vacation so I forgot to follow up. He needed an example that used solely SOFA code and I didn't have a chance to make a demonstration program that wasn't using astropy.

This came up again today in LSST where two times that were equal but went through UTC->TAI->UTC conversion and were no longer equal. Given that astropy knows that we have this picosecond error, is it possible that we can include this known error in the __eq__ implementation? Could we make the .precision attribute influence the accuracy of time comparisons? We are using nanosecond precision and don't care about a 10 ps error. This problem is making some people in LSST/Rubin question whether we should be using astropy at all (and instead go back to what we were doing which was integers with nanosecond precision).

@timj

@timj - I sent SOFA a patch years ago that would fix this round-trip error. I'll try to find it again and resend...

On __eq__, I don't think it's a good idea to change it; Time is effectively a high-precision float and I think it is fine for it to behave as such. But that also suggests a solution: I think it makes a lot of sense to have something like np.isclose for times; indeed, with the new __array_function__, it is quite easy to ensure that np.isclose itself can work on Time istances (see #8610; added this to its to-do list). Of course, for now you could do as we do in many of the time tests, np.abs(t1-t2) < ??*u.ns - explicit is better than implicit!

I'm 👍 on a Time isclose method, with a default that corresponds to our known expected precision.

Question is whether it is really like numpy np.isclose(a, b, atol, rtol) => astropy.time.isclose(a, b, atol) (where I suppose atol would be a quantity convertable to days), or a Time method Time.isclose(self, other, atol). I would generally prefer the latter as being more convenient (not requiring an extra import) but don't feel strongly.

I think the choice would be between np.isclose itself (via an __array_function__ override) or a method on Time. So, no extra import either way (assuming one has already imported numpy). Advantage of the former is re-use, of the latter having, e.g., a custom docstring.

p.s. I'd like the former to work regardless, as the meaning is obvious, with atol either a quantity with time units or a TimeDelta, and we can have the default be sensible. Extension to SkyCoord may be possible as well, since representations can be subtracted (but as yet do not have an __abs__ method - could be added, though; or just use .separation)

I think the choice would be between np.isclose itself (via an array_function override) or a method on Time. So, no extra import either way (assuming one has already imported numpy). Advantage of the former is re-use, of the latter having, e.g., a custom docstring.

There isn't any solutions for a custom docstring for the override case?

There isn't any solutions for a custom docstring for the override case?

Would be wonderful, but also quite tricky: how would one change the documentation on np.isclose? It does perhaps suggest a use for something like an astropy.time.np namespace that has documentation overrides. At a minimum, we'd need to document the numpy functions that are known to work (would be good for Quantity too...). It should be possible to auto-generate this, since there has to be a function that does the conversion to array before feeding it into numpy anyway (like https://github.com/astropy/astropy/blob/master/astropy/units/quantity_helper/function_helpers.py)

I do think users would quickly get used to the idea of just trying a numpy function that should make sense.

the shared common things could be inherited from the overriding __array_function__ itself? I really don't know how to resolve this though, but suspect it's a shared problem for all the libraries that are doing this.

Or maybe we just need it to be covered in the narrative docs, and as you say users will get used to it.

Looking at np.isclose, I'm a bit dubious of shoe-horning that to work for Time since rtol and equal_nan are both not applicable in this case. I would instead suggest to have the __array_function__ override raise NotImplemented and tell the user to use Time.isclose, which will have the proper documentation etc.

Sent an e-mail to Pat Wallace about the round-trip issue.

A new version of SOFA with improved precision is being prepared.

Good news. Of course, still really should have isclose...

New ERFA: liberfa/erfa#71

I'll take a crack on Time.isclose.

See #10646.

I think this is no longer an issue. This is what I get with astropy v5.2.1

>>> t3 == t1
True
>>> t2 < t1
False
>>> t1.jd2
-0.054203622508796234
>>> t2.jd2
-0.054203622508796234
>>> t3.jd2
-0.054203622508796234

Feel free to re-open if I misunderstood. Thanks!

Agreed that with .isclose() and the substantial work in ensuring rounding errors usually go the same way this is not really an issue any more. Thanks!

Navigation Menu

Repository navigation

Time equality comparisons are possibly too precise #8742

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Time equality comparisons are possibly too precise #8742

Description

Activity

mhvk commented on May 21, 2019

mhvk commented on May 21, 2019

mhvk commented on May 21, 2019

bsipocz commented on May 21, 2019

mhvk commented on May 21, 2019

bsipocz commented on May 21, 2019

mhvk commented on May 22, 2019

timj commented on May 24, 2019

timj commented on Jun 10, 2019

timj commented on Apr 7, 2020

mhvk commented on Apr 7, 2020

taldcroft commented on Apr 7, 2020

mhvk commented on Apr 7, 2020

bsipocz commented on Apr 7, 2020

mhvk commented on Apr 8, 2020

bsipocz commented on Apr 8, 2020

taldcroft commented on Apr 8, 2020

mhvk commented on Apr 8, 2020

timj commented on Jun 15, 2020

mhvk commented on Jun 15, 2020

timj commented on Aug 12, 2020

taldcroft commented on Aug 13, 2020

taldcroft commented on Aug 13, 2020

pllim commented on Jan 25, 2023

mhvk commented on Jan 25, 2023

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions