Time("J2000") and equivalents should be equal

I noticed this behaviour:

In [21]: t0 = time.Time("J2000", scale="tdb").replicate("iso"); t0
Out[21]: <Time object: scale='tdb' format='iso' value=2000-01-01 12:00:00.000>

In [22]: t0.jd1, t0.jd2
Out[22]: (2400000.5, 51544.5)

In [23]: t0.unix
Out[23]: 946727935.8160995

In [25]: t1 = time.Time("2000-01-01 12:00:00", scale="tdb").replicate("iso")
Out[25]: <Time object: scale='tdb' format='iso' value=2000-01-01 12:00:00.000>

In [26]: t1.jd1, t1.jd2
Out[26]: (2451544.5, 0.5)

In [27]: t1.unix
Out[27]: 946727935.8160994

In [29]: t1 - t0
Out[29]: <TimeDelta object: scale='tdb' format='jd' value=0.0>

In [30]: _.sec
Out[30]: 0.0

In [31]: t1.unix - t0.unix
Out[31]: -1.1920928955078125e-07

This is giving absolute differences of 1e-7. This difference may or may not be important depending on the specific application, but in any case are they expected? If yes, should they be documented? If not, what is going wrong?

Edit: I added the output of t1.unix - t0.unix to the summary.

WIP: Add propagate to date, test failing

@taldcroft

I'll let @taldcroft and @mhvk give a full answer but just to speculate on what might be happening here: t1.unix and t2.unix are just 64-bit floats, and -1.1920928955078125e-07 is basically the precision of the 64-bit floats:

In [2]: np.spacing(946727935.8160994)
Out[2]: 1.1920928955078125e-07

As such I don't know if the difference is significant. As you can see from:

In [29]: t1 - t0
Out[29]: <TimeDelta object: scale='tdb' format='jd' value=0.0>

the times are equal and the residual value for t1.unix - t2.unix is just a product of the inevitable loss in precision when going to 64-bit floats? This could also be related to the fact that in the original scale, the floats are exact (multiples of 0.5) whereas in the unix representation this is not the case?

@astrofrog

Thank you very much @astrofrog, makes total sense. However, this is the issue that originated this report:

In [26]: t0 = time.Time("J2000", scale="tdb").replicate("iso"); t0
Out[26]: <Time object: scale='tdb' format='iso' value=2000-01-01 12:00:00.000>

In [27]: t1 = time.Time("2000-01-01 12:00:00", scale="tdb"); t1
Out[27]: <Time object: scale='tdb' format='iso' value=2000-01-01 12:00:00.000>

In [28]: t0 == t1
Out[28]: True

In [29]: (t0 + 40 * u.min) == (t1 + 40 * u.min)
Out[29]: False

In [31]: (t0 + 40 * u.min)._time_difference(t1 + 40 * u.min)
Out[31]: 3.2337466038256935e-12

So perhaps a different algorithm in __eq__ should be used for comparison based on relative differences? Not sure how to do such a thing though, or if it makes sense at all.

It happens even with tiny additions:

In [45]: t0 == t1
Out[45]: True

In [46]: (t0 + 1 * u.nanosecond) == (t1 + 1 * u.nanosecond)
Out[46]: False

In [47]: (t0 + 1 * u.nanosecond)._time_difference(t1 + 1 * u.nanosecond)
Out[47]: -1.1546319456101628e-14

@taldcroft

I think this could be because internally one of the jd2 is much larger than the other, so then internally the two times have different precision. Actually looking at:

In [21]: t0 = time.Time("J2000", scale="tdb").replicate("iso"); t0
Out[21]: <Time object: scale='tdb' format='iso' value=2000-01-01 12:00:00.000>

In [22]: t0.jd1, t0.jd2
Out[22]: (2400000.5, 51544.5)

this doesn't seem right because having the two floats isn't improving the precision (they have similar orders of magnitude). @taldcroft - is that normal?

Hmm, the problem seems to be that the times are not being properly normalized to day and fraction upon initialization, and then by adding a number that has a small rounding error when represented as a fraction of a day, one gets a different result.

I think I see the problem: day_frac is not used in the TimeEpochDateString formats.

(I cannot make a fix right now - teaching in a few hours - but this should be quite trivial... In https://github.com/astropy/astropy/blob/master/astropy/time/formats.py#L1074,
write instead

jd1, jd2 = epoch_to_jd(iterator.operands[-1])
self.jd1, self.jd2 = day_frac(jd1, jd2)

Checking a bit further, I think the same needs to be done for TimeDatetime (just not to count on erfa doing it identically the same), and TimeString.

Since the fix is basically given in #6638 (comment) , I am tagging this for novice. Thanks!

@pllim

Thanks @pllim, @mhvk. I might fix this myself this week if I have the time.

@mhvk

@mhvk - hopefully I'm not misunderstanding, but it doesn't seem like the post-facto fix on jd1/jd2 will maintain full precision since jd2 is typically in the range of 50000 or 60000 instead of 1. So bits of precision are being permanently lost within the epoch_to_jd function (which is just a reference to the appropriate ERFA function).

Those ERFA functions are actually really trivial, e.g. the J2000 conversion C-code is just:

   *djm0 = ERFA_DJM0;  # Always 2400000.5
   *djm  = ERFA_DJM00 + (epj - 2000.0) * 365.25;

So maybe just rewrite this in Python in a way that guarantees jd1/jd2 are within the day_frac protocol?

I think the assumption in erfa is that the epochs are never that high precision. I think this is actually fine; I don't think there is any case where we would pass in a 2-element epoch. The problem with the sum being different for different input is that Time.__add__ assumes that jd1, jd2 are in day, fraction format, when they are not, and therefore the addition gives a slightly different result.

If we wanted, we could give a warning about possible precision loss when val2 is passed to any of the epoch routines...

@mhvk

@mhvk - fair enough. Maybe the easy fix will do the trick.

Accomodate jd1 and jd2, see astropy#6638

WIP: Sum test fails, see astropy#6638

Attempt to fix: #6653

WIP: Add propagate to date, test failing

Mark test as known failure

Navigation Menu

Repository navigation

Time("J2000") and equivalents should be equal #6638

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Time("J2000") and equivalents should be equal #6638

Description

Activity

astrojuanlu commented on Oct 1, 2017

astrofrog commented on Oct 1, 2017

astrojuanlu commented on Oct 1, 2017

astrojuanlu commented on Oct 1, 2017

astrofrog commented on Oct 1, 2017

mhvk commented on Oct 2, 2017

mhvk commented on Oct 2, 2017

pllim commented on Oct 2, 2017

astrojuanlu commented on Oct 2, 2017

taldcroft commented on Oct 2, 2017

mhvk commented on Oct 2, 2017

taldcroft commented on Oct 2, 2017

astrojuanlu commented on Oct 3, 2017

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions