Repository navigation
scalar Time + multidimensional Quantity = error #3039
Description
Adding a Time and Quantity works when the Quantity is a scalar and a 1D vector. If the Quantity array is reshaped to a 2D array, an OperandTypeError is thrown.
import astropy.coordinates as c
import astropy.time as t
import astropy.units as u
import numpy as np
base_date = t.Time('2014-01-01')
offsets = np.arange(0,10) * u.day
offsets
Out[7]: <Quantity [ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.] d>
base_date + offsets
Out[8]:
<Time object: scale='utc' format='iso' value=['2014-01-01 00:00:00.000' '2014-01-02 00:00:00.000'
'2014-01-03 00:00:00.000' '2014-01-04 00:00:00.000'
'2014-01-05 00:00:00.000' '2014-01-06 00:00:00.000'
'2014-01-07 00:00:00.000' '2014-01-08 00:00:00.000'
'2014-01-09 00:00:00.000' '2014-01-10 00:00:00.000']>
offsets.reshape( (2,5))
Out[9]:
<Quantity [[ 0., 1., 2., 3., 4.],
[ 5., 6., 7., 8., 9.]] d>
base_date + offsets.reshape((2,5))
---------------------------------------------------------------------------
OperandTypeError Traceback (most recent call last)
<ipython-input-10-4402dbceaf11> in <module>()
----> 1 base_date + offsets.reshape((2,5))
C:\Users\bnordgren\AppData\Local\Enthought\Canopy\User\lib\site-packages\astropy\time\core.pyc in __add__(self, other)
993 other = TimeDelta(other)
994 except:
--> 995 raise OperandTypeError(self, other, '+')
996
997 # Tdelta + something is dealt with in TimeDelta, so we have
OperandTypeError: Unsupported operand type(s) for +: 'Time' and 'Quantity'
Activity
cheekujodhpur commentedon Oct 19, 2014
Not too hard to solve by modifying the init for class TimeDelta. There are however two ways where your ideal output could go:
No.1 seems more intuitive to me. But I feel a discussion would help
bnordgren commentedon Oct 19, 2014
I was expecting No. 2 because (as a user), I thought that Time should be able to do anything that Quantity (and other numpy ndarrays) can do. In any case, the above error message is not at all helpful. I was tracking all over my code trying to figure out what broke this equation (which I didn't touch.)
The big win with No. 2 is the ability to write array equations which are the same regardless of dimensionality. For instance, a function or method which calculates something for a point does not need to be rewritten for use on a vector, and also does not need to be rewritten for use on a grid. For me, at least, I test equations using a point, then apply them to vector or grid data. Writing special cases for different dimensionality is irritating busywork, and leads to less maintainable code should you discover an error after you've implemented 3 versions. :)
My workaround is to reshape the offset array into a vector, then three or four steps later reshape a Quantity vector back into the original shape. This is OK if it can all be handled inside the same function/method, but is likely to lead to buggy code if the function/method has to return a Time vector without an indicator of the original shape.
The application, in this case, is an implementation of Meeus's rise/set/transit algorithm from "Astronomical Algorithms", 2nd ed. Each observing location (row) has three times (columns) associated with it. I discovered this issue when I tried to make the code work with a vector of observing locations instead of a single point.
mhvk commentedon Oct 19, 2014
@bnordgren - the problem here is that
Timedoes not allow multidimensional data:In principle, this should not be that difficult to change, as the times are internally stored as two-part julian dates,
jd1andjd2, which are bothndarrayand thus can have arbitrary dimension. The main issue is that it needs some thought on how to pass higher-dimensional arrays to theERFAroutines that do the actual conversion between timescales, etc. Still, with suitable.ravel()and.reshapecalls this should be do-able...@eteq, @taldcroft - will #2992 help with this?
cheekujodhpur commentedon Oct 19, 2014
Yes, I did not realize the suitability to specific applications.
I would like to start workin' on this, testin' a few ideas while we wait for more input
mhvk commentedon Oct 19, 2014
@cheekujodhpur - would be great! The actual check is done in
_make_1d_arrayinastropy/time/core.py. Unless the newerfaroutines can take multidimensional arrays (perhaps unlikely), all the calls toerfaneed to be wrapped. As an example, the first place would be athttps://github.com/astropy/astropy/blob/master/astropy/time/core.py#L389 - where one would probably replace the line with
and then before the new time is initialized, at https://github.com/astropy/astropy/blob/master/astropy/time/core.py#L409, set the old shape back, as in:
and similar changes in other locations where
erfais usedtaldcroft commentedon Oct 20, 2014
In principle I'm fine with relaxing the restriction on having 1-d arrays. The original motivation was that it kept things simple and that we could wait until someone actually raised the issue. If there is a need, and someone to do the implementation then it sounds good.
eteq commentedon Oct 20, 2014
@mhvk - I don't think #2992 will impact this too much - this is more about how
Timeviews itself, and ideally #2992 will be a transparent change fortime.mhvk commentedon Oct 20, 2014
@eteq - but I see discussions about broadcasting in #2992, which, if this is true broadcasting as for normal numpy operations, implies there will have to be no changes to
Time. Am I misunderstanding?mhvk commentedon Oct 21, 2014
@cheekujodhpur - now looked at #2992 a bit more and I do think this will solve the broadcasting for us, so I would hold off on experimenting too much here.
mhvk commentedon Nov 21, 2014
@cheekujodhpur - With #2992 merged, I think it will now be quite trivial to ensure
Timecan handle multidimensional data; in fact, all that may be needed is to remove the check on the dimension! Let me know if you are you still interested in trying this (if not, I'll try myself).mhvk commentedon Nov 24, 2014
@cheekujodhpur - I started with converting
Timeso that scalars are dealt with better and it turned out that this also made the multi-dimensions work. Should have a PR with it shortly, but thought I would let you know so that we don't duplicate each other.mhvk commentedon Nov 24, 2014
p.s. See #3135 for an ERFA bug that needs to be fixed before I can send my PR...
mhvk commentedon Dec 18, 2014
This now is solved by #3138, so I'll close this. @bnordgren, thanks for reporting this issue, it will no longer be one in 1.0! @cheekujodhpur, thanks for your offer to help, even if here it ended up being taken over by events.