Add transform method to all representations

[nstarman]

Description (edited)

All representations now have a transform method, which allows them to be
transformed by a 3x3 matrix in a Cartesian basis. By default, transformations
are routed through CartesianRepresentation. SphericalRepresentation and
PhysicssphericalRepresentation override this for speed and to prevent NaN
leakage from the distance to the angular components.
Also, add a utility function in matrix_utities to check whether a matrix is
a rotation (proper or improper).

@mhvk @adrn
closes #11423
Blocked by #11568

Signed-off-by: Nathaniel Starkman (@nstarman)

[pep8speaks]

Hello @nstarman 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

There are no PEP8 style issues with this pull request - thanks! 🎉

Comment last updated at 2021-04-20 17:16:15 UTC

[nstarman]

Now I'm looking at BaseAffineTransform to see how to route transformations through this new machinery. I'm thinking of re-implementing the rotation matrix check if the representation is SphericalRepresentation. If it is, that transform method is called over the Cartesian method. To prevent repeating the rotation matrix check, SphericalRepresentation.transform accepts a flag, so that can jut be passed by BaseAffineTransform.

[mhvk]

General comment: I don't think we use a non-cartesian matrix anywhere, and I must admit I have a hard time seeing what the use would be - in fact such a matrix cannot currently be represented properly as a Quantity. I would suggest focussing on just implementing the cartesian variant here, which does have direct use.

[mhvk]

My sense would be to use this PR only to change the representations and test those. In particular, I'd make a general transform method on BaseRepresentation, which just goes through Cartesian (i.e., does what is now done inside the transformation mechanism) and then override it on Spherical and UnitSpherical (have to do the latter too, since the general case will no longer have unit distance).

In the transformation machinery, in the first instance we can then simply stop worrying about converting to/from cartesian, and in the second instance we can start to call using the explicit flag for cases where we know a priori it is a rotation matrix.

[mhvk]

On the spherical matrix for transform, maybe another way to express why I'm not quite sure that is useful is just to think that in the end one if multiplying the longitude and latitude with numbers, which does not seem to be very useful!

[nstarman]

@mhvk, I reimplemented the representation transformations as you suggested.
It works! the only thing that doesn't work is if the matrix is not (3,3), but (N, 3, 3) like in some of the tests when a parameter like "obstime" makes the matrix broadcast. I need to fix the unitary matrix check for that edge case.
It even works for Coordinate Frames, but it's a bit inefficient because each called to transform the representation performs a unitary matrix check (because the flag "is_rotation" is None). I'm working on propagating the flag up the coordinate transform machinery so this can be set for known rotations.

[nstarman]

I've implemented the flag "is_rotation" for static and dynamic matrix transforms.

[nstarman]

@mhvk @adrn , rotation matrix checks now broadcast. All built-in coordinate transforms that use a StaticMatrixTransform and DynamicMatrixTransform now use the flag "is_rotation". #FightTheNaN

Edit: I'm having trouble diagnosing the errors. Many seem related to AffineTransform, which should not be impacted by this PR.

[mhvk]

@nstarman - while reviewing your new additions, it suddenly struck me that perhaps we are making it too complicated. How about we we do define a general .transform but instead of trying to track whether something is a rotation matrix, for Spherical we simply split off the distance piece, i.e., we do the equivalent of

def transform(self, matrix):
    return self.represent_as(UnitSpherical).transform(matrix) * self.distance

(in principle, to avoid this being slower, one might implement this a bit more explicitly, e.g., using erfa.s2c as in UnitSphericalRepresentation.to_cartesian()).

With that, the only other thing needed would be to remove the transformations to Cartesian in the coordinate machinery. What do you think?

[mhvk]

Hmm, had some comments, but am now thinking we should just write the transformation such that the distance is done separately anyway - see in main thread... (so, this is not quite complete either...)

[mhvk]

@nstarman - before we go further, what about my larger suggestion, at #11444 (comment)? It avoids even the need of checking/tracking whether something is a rotation matrix.

[nstarman]

@mhvk, I really like the idea, and it works wonderfully for representations, but I don't think it works for differentials.
In my attempts to get this to work I found that while unit representations can be rescaled multiplicatively, this is not true for the differentials. The unit differentials has a radial component of 0, not 1, so the "rescaling" must be additive.

Consequentially, if the transformation matrix is not a rotation, setting the radial component to 0 before applying the matrix destroys information. Therefore, only when the radial effect of the matrix is known a priori (eg none for a rotation) can we shortcut in this manner. This can be seen in the following code snippet, which is the transformation for differentials in the SphericalRepresentation.transform, when the matrix is a rotation matrix: for the representation we can multiply by the distance (or assign since it's 1), for the differentials we must add (or assign since it's 0).

Edit: as a more detailed demonstration, the radial component of the differential cannot be rescaled and the process is therefore not information-preserving.

I haven't looked into speedups using erfa. Perhaps going directly to the matrix manipulations would circumvent this issue...

[mhvk]

@nstarman - I think it could still work in principle, if one took care to split out the radial velocity component and treat that separately (for higher order derivatives this becomes quite painful, but we cannot deal with those properly yet anyway).

But perhaps it still makes sense to decouple the representations from the differentials? First move .transform to all the representations, and ensure things still work when the change to/from cartesian is done there, then implement the distance-avoiding change for Spherical, just for the case of no differentials. And only after that worry about whether there are other corner cases that can be dealt with (and whether those need a is_rotation flag). This PR could then avoid dealing with differentials altogether, keeping things simple (and if we avoid 4.3, we can ensure that we do not bake in any decision on extra keyword arguments). And of course the main benefit is that we can go in small steps...

[nstarman]

Note: #11470 will simplify rescaling. Don't need to strip differentials.

[adrn]

Overall looks great. Mostly minor nitpicks or questions from me!

[mhvk]

So, with #11568 the PhysicsSpherical tests should all work too. Shall we wait for that?

[nstarman]

🤞

[mhvk]

Looks all OK, but on final review one super nitpicky comment about the tests... If you're fed up, though, I'll merge as is...

[adrn]

Would ideally like to take another look at this, but since you addressed my first round of comments, @mhvk has taken another close look, and I want to see this in...approving!

[mhvk]

And let's get this in as well! Next, the differentials... (but maybe after #11470?)

[nstarman]

Wooh! Thanks for the collab @mhvk @adrn. Good to get this in.
Agreed on the differentials and for waiting until after #11470.

It seems the overall plan of action is: #11470 -> .transform on the differentials -> shortcuts including differentials.
A good complementary objective is minimizing intermediate representations in CoordinateFrame objects.

[nstarman]

Should this be in the whatsnew?

[mhvk]

My sense is to have a what's-new together with the larger overhaul in #11470 - which means perhaps for 5.0 - since the whole makes working with representations as vectors much easier. On its own, there is probably no reason for the average user to know about this, as it will mainly have (nice!) impacts on our internal implementation of coordinate transformations.