Allow TimeDelta to have (well-defined) scales other than TAI - #1932

@taldcroft - above a first attempt; one particular item I'm unsure about is what to do with conversion to/from Quantity. As is, I completely ignore the scale, which is consistent with current initialisation, and perhaps most logical given that if one has a TCB, one presumably wants to work in TCB seconds as well. Still, it may lead to confusion.

Reworked it a bit, trying to minimize scale changes, and adding UT1 as a possible TimeDelta class that is cannot be converted to other classes (but is well-defined enough on its own, though maybe should be excluded for conversion to/from quantity).

@mhvk - overall this looks really good. The main thing I'm struggling with now is that it seems intricate. I'm just thinking ahead to documentation and whether there will be surprises for users. There is also a component of code reliability, but I'm sure this can handled with testing.

I've been trying to think about ways to simplify the behavior (and code) by imposing some restrictions. To be honest I haven't converged to anything that is really great, but I'm still thinking about it. An extreme example is requiring the same scale for any arithmetic. At some level this makes everything more explicit and removes a layer of magic that might lead to problems. On the other hand it takes away some of the nice "it just works" feel. Hmmm.

@taldcroft - I'll discuss both UT1 and the general case here, so we don't "loose" it if files get updated.

I've also been struggling with what to do, and my reasoning ended up being based mostly on whether a time difference makes sense if one only knows the scale, i.e., I tried to take the reasoning that led to geocentric and barycentric groups to its logical conclusion.

With this condition, UTC has to be excluded because for a given dt_utc=utc2-utc1 one cannot guarantee that ut2==ut1+dt_utc without storing ut1 and ut2 or some other extra information in dt_utc. In contrast, for all other timescales, one can make this guarantee within the timescale itself, and it requires only a well-defined transformation if one stays within the three groups of timescales defined here. Hence, I thought one might as well include UT1.

Of course, UT1 also has direct physical meaning: even if the clock doesn't quite run at a constant rate, dt_ut1 is still a well-defined time representation of a rotational phase (which of course time was for millenia). In this respect, the difference with TAI seems similar to the the difference between TAI and TDB (for one you need to know as-yet-incomplete tidal theory, for the other GR).

On your other question of whether one should just stay within a given timescale: to me, UTC kills this; I feel one should be able to take time differences in UTC, since that will be the most common case, but for those the time scale has to change to TAI. Given that, we might as well allow all changes within a group.

In the end, perhaps the right question is whether the behaviour will generally be as expected by the user, and whether, what is not, is easy to explain in the documentation. I think it is OK: for most timescales, say, a second is a fixed fraction of the Cs hyperfine line frequency, but it depends on location, for UT1 it is 1/86400 of a solar day.

Maybe the next thing I should do is try to adjust the documentation and add further tests?

What if I just wander down a path and see where I stumble...

Imagine (as I mentioned) that TimeDelta is a purely mathematical construct that simply represents the difference in JD of two like-scaled times. It would seem that dt_utc = utc2 - utc1 is guaranteed to be transitive with utc2 = utc1 + dt_utc since you are just manipulating the JD values (which is a continuous variable in the expected way even for UTC).

To me there is no difference between UTC and UT1 in this respect. Apart from the heritage of UT1 as angle, they are both just timescales where the physical duration of a day or second varies with time. For UTC it's a step function on one day every so often, for UT1 it's continuously variable. In both cases I think it's natural that if you add 1.0 JD in that scale then you come back to exactly the same 24-hour wall time. That would be true in this scenario. I think this is an easily explainable and consistent definition, and it's very nearly what you've already done.

The only potentially unexpected consequence would be, e.g. adding 0.5 JD on a day with a leap second, in which case you accrue a fractional bit of second in the ISO representation.

I think that one nice consequence of this is that TimeDelta(xx, scale=None) could be allowed and would make sense as the universal time delta which would be allowed for arithmetic with a Time in any scale. It would be exactly the same as a time-like Quantity.

Part II next...

Here what would makes sense to me for the rules of combination, without reference to what's done so far (but I think it mostly matches). I'm using the nomenclature T(X) as a Time object with scale X, and dT as a TimeDelta. In some cases I say "add" when I mean "add/subtract".

• T(X) +/- dT(X)
  Add JDs.
• T(X) +/- dT(None) or T(X) +/- Quantity(time-like)
  This is always allowed and just adds JDs. A day is defined as 86400 seconds.
• T(X) +/- dT(Y)
  Do (in essence) (T(X).Y +/- dT(Y)).X. In practice the dT value could instead be converted to scale X if possible. By doing this, if T is convertible to scale Y then the operation always succeeds. Return value always has scale X.
• T(X) - T2(Y)
  Do T(X) - T2(Y).X and return dT(X).
• dT(X) +/- dT2(Y)
  Add JDs for dT(X) and dT(Y).X. I.e. this will succeed if dT(Y) can be converted to scale X via the SCALE_OFFSET factors. Returns delta time in scale X.

Addendum:

It might be good to disallow adding the variable time scales UTC and UT1 to anything but themselves.

@taldcroft - nice idea to have TimeDelta(xx, scale=None) as a more precise equivalent to using a Quantity with time units! In my implementation, I felt it was rather awkward that, effectively, time-like quantities had arbitrary scale, but with this, it becomes part of a larger scheme.

I think your logic for what should happen matches mine, except that I do not think we can handle UTC this way. The example with 0.5 JD is a good one: essentially it is not possible to always have that 1 day == 86400 second. The broader problem is perhaps easiest seen for T(UTC) - T2(UTC):

In [8]: Time('2012-06-29T00:00:01', scale='utc').jd2 - Time('2012-06-29T00:00:00', scale='utc').jd2
Out[8]: 1.1574074074074073e-05

In [9]: Time('2012-06-30T00:00:01', scale='utc').jd2 - Time('2012-06-30T00:00:00', scale='utc').jd2
Out[9]: 1.1573940116433837e-05

So, in your scheme, we would get different answers just by being on the day of a leap second. I think very few users would expect this, so see no alternative to sticking with your original scheme of converting to TAI for this case (I just cannot bear to think the number of wrong periods, etc., that would be published!). Note that for UT1 this problem does not occur: the second may differ in atomic time for different epochs, but there are always 86400 of them in a solar day. I also do not worry as much about user mistakes, since direct usage of UT1 is rare and already implies some awareness of timescales.

Note that for UT1 this problem does not occur: the second may differ in atomic time for different epochs, but there are always 86400 of them in a solar day.

In this scheme there are also exactly 86400 UTC seconds in each UTC day, by definition. It just so happens that on leap second days the UTC second is a bit longer (86401 / 86400 TAI seconds). The math all works, and essentially we accept this already because in ERFA, 12:00:00.000 hours is not 0.5 JD on some days.

Having said this, I see your point that this could very easily trip up users that aren't paying attention. I think the same could be said of UT1. What about this:

• Any time delta T1 - T2 where T1 is UTC or UT1 returns the TimeDelta in TAI.
• Remove the UT1 TimeDelta scale since it is a potential source of errors. If users complain then reconsider our options.
• Define time-like Quantities to be TAI, always. This gives a single unambiguous definition of what that means which is most likely to be what people expect. People that really care about time wouldn't be using Quantity anyway and with this PR they can manipulate Times conveniently within their scale of choice using TimeDelta.
• Require a scale for TimeDelta object creation. Plan-B would be having the default be TAI by analogy with Quantity.

So, in your scheme, we would get different answers just by being on the day of a leap second.

The "different answers" bit depends on your representation of time (aka the format). If you are using JD as the format, then adding 1 "second" (by which you meant incrementing the "second" value in the ISO representation, i.e. a TAI second here) will result in different JD answers if you are on a day with a leap second.

@mhvk , @taldcroft - Do you think you have reached any agreement on what to do for this? If you want my two cents:

I don't think there is any problem with allowing UT1 TimeDeltas. In my opinion, if someone is using UT1 times they are likely doing so on purpose, and would expect t2.ut1 - t1.ut1 to a) be possible, and b) return a value in UT1 seconds, not in TAI or some other scale. Essentially this is the same as my experience trying to get TDB differences. It's hard for me to imagine a case where someone unwittingly ends up working in UT1 and runs into errors because of this.

UTC is a bit trickier, because I can definitely imagine a more casual user working with UTC times who wants differences in TAI, without even necessarily realizing there is an issue... so I think there is a valid argument for auto-converting UTC differences to TAI. On the other hand, the "purely mathematical" scheme proposed above by @taldcroft where time differences are always just the numerical difference of the jd values does have a nice logical consistency to it. I could go either way on this one.

Finally, aside from a possible special case for UTC/TAI, auto-converting between scales (ie allowing "T(X) +/- dT(Y)" and similar) still seems like a potential source of confusion/bugs to me. I think the simplest thing to do would be to raise an exception in cases like this, which should force the user to explicitly state which scale to use for the calculation.

@demorest - I agree we should move forward on this. I've actually been thinking about it, and pondering putting in a PR on @mhvk 's branch to at least demonstrate what I was proposing. It can be easier to see how something plays out with code in hand.

I don't think there is any problem with allowing UT1 TimeDeltas. In my opinion, if someone is using UT1 times they are likely doing so on purpose, and would expect t2.ut1 - t1.ut1 to a) be possible, and b) return a value in UT1 seconds, not in TAI or some other scale.

What about allowing subtraction to form a TimeDelta in UT1, but not allowing any subsequent arithmetic operations with that TimeDelta instance? To put it another way, what is the scientific use case for adding a UT1 Time and a UT1 TimeDelta, when it is known that this will give a meaningless answer? Would any of the expert users that are using UT1 ever need to do this, and why?

@taldcroft - I have no idea why someone would need to do a lot of UT1 arithmetic. But as I think has been discussed already in this thread, I don't see why adding a UT1 TimeDelta to a UT1 Time is meaningless. It seems to me that it's a well-defined operation. So why not allow it?

Finally, aside from a possible special case for UTC/TAI, auto-converting between scales (ie allowing "T(X) +/- dT(Y)" and similar) still seems like a potential source of confusion/bugs to me. I think the simplest thing to do would be to raise an exception in cases like this, which should force the user to explicitly state which scale to use for the calculation.

@demorest - In one sense I'm with you here to be strict and conservative, and it's easy enough to raise an exception with a helpful warning that says (effectively) "You've tried to do arithmetic with two objects that have different time scales. In order to do this, you need to convert them to the same scale." But, at the same time I often get annoyed at libraries that behave like this, because astropy.time knows exactly what the problem is and it knows how to fix the problem. So why not just do it?

The question is then: what might go wrong? It seems to me that if we have the appropriate metadata to get T(X) over to scale Y and back to X then everything is fine. Am I being naive or missing something subtle, perhaps with the GR scales?

Here is my understanding, please correct if I have something wrong. What I'm saying is likely blindingly obvious, I'm just trying to write things out in detail.

UT1 is unique among the astropy.time supported time scales. For all the others, in the local frame of reference corresponding to that scale, a duration of one second in that scale corresponds to a fixed number of transitions of an atomic clock. This is independent of the epoch, so that a TimeDelta of one second is always a fixed physical time in that reference frame. For any time T0, if you compute T1 = T0 + 1.0 sec, the number of atomic ticks between T1 and T0 will be constant. (Again, the observer and clocks are all within the relevant local frame of reference).

Instead if you use UT1 to determine the orbital period P of planetary moon from 1970 data, and then you predict the exact transit time in 2014 by computing something like T1 = T0 + P, you get the wrong answer. Using time for astronomical analysis generally implies physical time.

So disallowing arithmetic like T0 +/- dT (for both in UT1) is applying a small degree of conservatism by saying "Hey, what you are doing is likely wrong if you are trying to derive a physical result like the expected time of an event".

We can put in this restriction, and even include text in the exception instructing the user to mail astropy-dev if they really have a use case. It's always easier to add back features than to remove them later.

Well, my take on this is that there is nothing unphysical about UT1. It's just based on Earth rotations rather than frequency of atomic transitions. I don't think the moon transit example is a good argument for disallowing UT1 arithmetic. I could rewrite this example to say how doing millisecond pulsar arrival time analysis in TAI gives you the wrong answer (it does and this is why we use TDB or TCB). But that doesn't mean TAI should never be used, it just means different time scales have different purposes.

In practice of course the Earth is not a very good clock, and the main purpose of UT1 is to figure out where your observatory is at a given time, not to compute elapsed time between events. I really can't see how someone would accidentally end up doing their moon transit analysis (or whatever) in UT1... but I'm also not sure that it's astropy's job to stop someone from doing this if they wanted to. Disallowing UT1 math does put a roadblock in place for someone who has a legitimate reason to use this feature (again, I don't know specifically what this use would be but that's not really the point). To me this seems like a much more annoying example of "overprotective software" than not having automatic timescale conversions.

In the end, however, the main thing I care about is doing high-precision TDB differences. So whatever gets that implemented fastest will be fine with me ;) And again, thanks for all your work on this, it really is a very useful package!

I could rewrite this example to say how doing millisecond pulsar arrival time analysis in TAI gives you the wrong answer (it does and this is why we use TDB or TCB). But that doesn't mean TAI should never be used, it just means different time scales have different purposes.

Your point of regarding UT1 as being derived from a lousy clock, and it's purpose of computing observatory angles, is well taken. So I drop my objection to arithmetic with UT1 times. Thanks for your patience, I'm starting to understand time a little better (having previously understood only just enough in order to wrap SOFA).

@demorest, @taldcroft - great to see convergence, especially since it seems towards what I had done for this PR 😃 I think I still like @taldcroft's summary of what the difference additions/differences would do, and will work with that (including trying to support TimeDelta(..., scale=None)). @taldcroft - after you made the list, you suggested being less flexible; let me know if you still feel that should be the case.

@mhvk - when you said "being less flexible" was that referring to T(X) +/- dT(Y)? If so, unless someone shows a concrete example of how that could go wrong, I would keep that.

@taldcroft - yes, that, and to perhaps not have TimeDelta(..., scale=None). I think I'll start with the more flexible approach, as it is closer to what I have. As usual, it may also help to adjust the documentation -- often, I only find out that something is really quite illogical when describing it!

👍

Any news on this? It is impacting (in a minor way) decisions we are making in PINT design... Thanks!

@scottransom - sorry, have been distracted by other things. The decision is there, just needs implementation (likely not that different from what's here) and that needs me (or someone else) to find time.

No worries.... I was just kicking the tires a bit... :-)

@taldcroft - I made all the changes, mostly staying on the "allow as much as possible" side (following your logic; indeed, tried to explicitly test these in the new test_time_delta.py). Also made a first cut at additional documentation.

@taldcroft - I got travis to pass also the tests in the documentation now. What do you think? OK to merge?

@mhvk - this looks awesome. So the only final quibble is to run pep8 on core.py and fix the warnings. I use this PEP8 config:

neptune$ more ~/.config/pep8 
[pep8]
max-line-length=100
select=E,W

neptune$ pep8 astropy/time/core.py
astropy/time/core.py:89:41: E226 missing whitespace around arithmetic operator
astropy/time/core.py:89:45: E226 missing whitespace around arithmetic operator
astropy/time/core.py:91:42: E226 missing whitespace around arithmetic operator
astropy/time/core.py:91:46: E226 missing whitespace around arithmetic operator
astropy/time/core.py:93:42: E226 missing whitespace around arithmetic operator
astropy/time/core.py:93:46: E226 missing whitespace around arithmetic operator
astropy/time/core.py:197:68: E226 missing whitespace around arithmetic operator
astropy/time/core.py:524:50: E226 missing whitespace around arithmetic operator
astropy/time/core.py:528:50: E226 missing whitespace around arithmetic operator
astropy/time/core.py:1451:12: E111 indentation is not a multiple of four
astropy/time/core.py:1453:12: E111 indentation is not a multiple of four
astropy/time/core.py:1455:35: E127 continuation line over-indented for visual indent
astropy/time/core.py:1945:61: E226 missing whitespace around arithmetic operator

Whitespace around operators makes the code easier for me to read.

@taldcroft - OK, made the PEP8 changes -- including ensuring line lengths are < 80 characters (my emacs window is that size, so this is something I actually like...). The only ones I skipped where
...wrap_angle=180*u.degree... in longitude initializers, where I think spaces around the * would make the code less readable.

OK, you're go to merge! Thanks @mhvk, this is a big improvement.

Yes, thanks @mhvk! I'll try to play with this in PINT on my long flight home today...