Repository navigation
Problems with basing CartesianPoints on plain 3xN Quantity arrays #1850
Description
While looking at the test cases for Coordinates to find out why changes to Quantity I made caused problems (#1848), I noticed this in test_distance.py:
#always preserve the CartesianPoint's units
c5 = 3 * u.pc + c
assert c5.unit == c.unit
which is really odd physically: what coordinate should the 3 pc be added to?
This brought home something I worried about while looking for the problems: by using a 3,N array, how can CartesianPoints protect against unphysical mixing of x, y, and z? Or ensure that methods such as np.mean work only work along relevant axes?
Easiest may be to behave more closely to the SphericalCoordinatesBase and just have x, y, and z Quantity arrays internally, but if they are combined in a single array, my suggestion would be to use a record array as the basis instead (with dtype=[('x', 'f8'), ('y', 'f8'), ('z', 'f8')]). In either case, one needs to override methods for the relevant functions.
Activity
embray commentedon Dec 2, 2013
On the question of
np.mean, one can already use theaxis=argument to specify which axis to take the mean along. Or if you just want the mean of one coordinate you can take the mean of a slice. But I agree it's not a bad idea to use a record array. Not sure how that works with Quantities but it should be okay in principle...mhvk commentedon Dec 2, 2013
@embray - indeed, the current implementation does not prevent one from doing the right thing (though one has to keep remembering the order in the array), but ideally it should prevent one from doing something as unphysical as the test case... Fortunately, it does seem
Quantityis quite happy with record arrays (though perhaps subclassing fromTableis not completely crazy either; that way spherical coordinates could have the same structure (three columns as well, but with different units).eteq commentedon Dec 12, 2013
@mhvk - there are three reasons for this that I remember off the top of my head:
x,y,zarrays and re-forming them every time, at least for some array sizes.cin the above example is a 3x1 array, it's the same as offsetting the origin by that vector. With a record array form, you'll always have to completely re-build the array instead of letting numpy optimize this.Quantitybeahvior just transfers over. Although I do recognize your point is that sometimes this is a bad thing if the behavior seems not as applicable. To me consistency makes things like this easier to understand, and therefor has a slight advantage over it being "exactly right". But that's all pretty vague, so there's definitely room for debate and/or change.One other thing: it may be that this should be part of the upcoming new coordinate API/APE (it's been a bit delayed because I decided I should do an APE on the topic at the same time). CartesianPoints will probably be something of a prototype for other similar classes there, so perhaps that context will help clear some of this up.
mhvk commentedon Dec 12, 2013
@eteq - Indeed, for performance reasons you're almost certainly right that separate arrays is a bad idea.
recarrayis much less bad, since one can quickly change it to a normal array just by doing the appropriateview. But your other points are good too, and probably it is indeed best to think of this more in the larger context. In this respect, one possible plus for therecarrayis that it may be substantially easier to share code withSpericalCoordinates(which might have an'r','phi','theta'recarray inside).eteq commentedon Dec 12, 2013
@mhvk - this is a good idea that definitely bears more thought... I'll have to do some performance tests to see how much (if any) of a penalty `recarr.view(dtype=float).reshape(...)`` has. Will definitely want to think more about it for the upcoming API!
astrofrog commentedon May 9, 2014
@mhvk @eteq - is this still relevant post-APE5?
mhvk commentedon May 9, 2014
No, it isn't; feel free to close (am offline now; sending a reply e-mail just so this doesn't get forgotten).
eteq commentedon May 9, 2014
Agreed, so I'll close this just so we don't forget to do so after merging APE5.