Representation arithmetic with differentials by mhvk · Pull Request #11470 · astropy/astropy

mhvk · 2021-04-03T01:41:38Z

~~Note: built on top of #11469, so that should be dealt with first! (Hence making this a draft, even though it is ready otherwise.)~~ (edit: #11469 is merged and this is rebased)

With this PR, differentials are taken along for multiplication, division, and negation:

In [3]: c = CartesianRepresentation([1.,2.,3]*u.kpc, differentials={'s': CartesianDiffer
   ...: ential([-2., 3., 4.]*u.km/u.s)})

In [4]: c2 = c * 10

In [5]: c2
Out[5]: 
<CartesianRepresentation (x, y, z) in kpc
    (10., 20., 30.)
 (has differentials w.r.t.: 's')>

In [6]: c2.differentials
Out[6]: 
{'s': <CartesianDifferential (d_x, d_y, d_z) in km / s
     (-20., 30., 40.)>}

For all other representations, the operations will continue to be equivalent to transforming to cartesian, applying the operation, and transforming back.

fixes #10987

p.s. No problem if this does not get in 4.3.

github-actions · 2021-04-03T01:42:19Z

👋 Thank you for your draft pull request! Do you know that you can use [ci skip] or [skip ci] in your commit messages to skip running continuous integration tests until you are ready?

mhvk · 2021-04-21T13:58:54Z

@nstarman - maybe could you have a look at this PR?

nstarman · 2021-04-21T18:11:53Z

        # super() checks that the class is identical so can this even happen?
        # (same class, different differentials ?)


It doesn't look like BaseRepresentationOrDifferential checks equality of differentials, so this check is necessary. Fixing the f-string does mean this line is counted as not being tested... Maybe add a quick test? Or out of scope.

nstarman · 2021-04-21T18:37:14Z


+        if self.differentials:
+            for key, differential in self.differentials.items():
+                diff_c = differential.represent_as(CartesianDifferential, base=self)


I was thinking about this for a while, but I think this is the correct route: don't assume that a subclass has defined _scale_operation and default to routing through Cartesian.
My one suggestion, and apologies because it's a bit of a bigger change, is to not do the scaling here, but delegate it to the method on the differential. Looking at CartesianDifferential, it inherits _scale_operation from BaseDifferential, which scales all components. However, ALL the other differentials also inherit the method... isn't this a bug? I would suggest moving the current _scale_operation from BaseDifferential to CartesianDifferential and implementing on BaseDifferential a similar method to here, where things are re-represented as a CartesianDifferental and then scaled appropriately. I think this would a) fix a bug and b) mean that as _scale_operation methods are added to each differential, they would automatically be used here.

I'm coming back to this only now. And remember that I struggled with this as well. In the end, it is, perhaps oddly, not a bug that _scale_operation just scales things for all differentials. E.g., consider

sph_coord + sph_diff * 1 * u.yr

In this case, the proper motions in sph_diff should just be multiplied by the time directly - there is no need to transform them first.

But of course, this is not correct if one just does sph_rep * 10 - then the angular components should not be affected. So, I went with your suggestion to move this (mostly) up to the differentials, but they need to know whether their base was scaled as well. It is all becoming a bit messy (maybe like your .transform) and perhaps there is a better solution...

That's a good point about unit multiplication. As you point out, there is a need to differentiate between types of scaling. What I'm currently leaning towards for .transform is a private argument to provide the scaled based.

Now I see the scaled_base argument. Agreed on the messiness of needing to track whether base was scaled. At least it's a private method and can change.

nstarman · 2021-04-21T18:40:43Z

-        self._raise_if_has_differentials('division')
-        return self._dimensional_representation(lon=self.lon, lat=self.lat,
-                                                distance=1. / other)
+    def _scale_operation(self, op, *args):


If the distance stays 1, I think it stay a UnitSpherical. Correct me if I'm wrong, but it doesn't appear to do that in SphericalRepresentation. I do see a #TODO, is this what that refers to?

It only gets here for multiplication and division; I felt it was a bit too much to check whether other is identical to 1.

nstarman · 2021-04-21T18:42:05Z

    def __neg__(self):
-        self._raise_if_has_differentials('negation')
-        return self.__class__(self.lon + 180. * u.deg, -self.lat, copy=False)
+        if any(differential.base_representation is not self.__class__


When does this happen?

A UnitSphericalRepresentation can carry a SphericalDifferential or a RadialDifferential. Quite messy...

But clearly not tested - now done.

I think a short comment would be good.

nstarman · 2021-04-21T18:44:08Z

+               for differential in self.differentials.values()):
+            return super().__neg__()
+
+        result = self.__class__(self.lon + 180. * u.deg, -self.lat, copy=False)


Probably for a followup, but do we want to try to stay on the same branch cut? One thing I've been toying around with, to little success so far, is trying to keep a memory of the phase wrap. Useful for dynamics and I'm sure other fields.

I think the class will wrap it, but, yes, let's do it as follow-up.

nstarman · 2021-04-21T18:57:39Z

-        return self.__class__(self.lon + 180. * u.deg, -self.lat, self.distance,
-                              copy=False)
+    def _scale_operation(self, op, *args):
+        # TODO: expand special-casing to UnitSpherical and RadialDifferential.


e.g. reduce back to Radial or UnitSpherical, if appropriate ?

This will reduce back to Radial or UnitSpherical, and, following your suggestion, no longer goes through Cartesian, so there is not much of an effect any more.

nstarman · 2021-04-21T19:05:05Z

Oops. Forgot to leave this comment when I pressed submit review...
I'm really looking forward to having this built in! My main concerns are related to _scale_operation on differentials and how this might avoid a bug, but also allow for delegation when scaling differentials in representations to the differentials themselves.

eteq · 2021-05-03T14:59:39Z

Since this hasn't made the feature freeze cutoff and is not a bug/docfix I'm re-milestoning to 5.0

nstarman · 2021-05-04T00:29:36Z

        except Exception:
            return NotImplemented

+        if self.differentials:


As I noticed when working on .transform stuff, isn't self.differentials always a dictionary, even if it's empty?
In which case...
for key, differential in self.differentials.items(): will just skip if there are no differentials. The loop doesn't need to be protected.

Good point! Fixed.

nstarman · 2021-06-01T18:15:23Z

Let me know how I can help.

mhvk · 2021-06-01T19:52:13Z

@nstarman and @adrn - I think for this PR the main question is whether the addition of scaled_base is fine, or whether there is a better solution. Since it is a private method, at some level it doesn't matter - we can change it later anyway.

@nstarman - one thing that may be helpful for @adrn is perhaps resolve conversations where you feel we've reached a decision.

nstarman · 2021-06-01T19:56:12Z

@nstarman - one thing that may be helpful for @adrn is perhaps resolve conversations where you feel we've reached a decision.

I actually don't have the GH privilege to resolve conversations except on my PRs. @pllim is there something to which I can apply for these privileges?

nstarman · 2021-06-01T19:59:56Z

@nstarman and @adrn - I think for this PR the main question is whether the addition of scaled_base is fine, or whether there is a better solution. Since it is a private method, at some level it doesn't matter - we can change it later anyway.

I think it's fine, largely because it is a private method. Obviously the transformed_base argument in Differential.transform() mirrored this PR, but since that is now merged, scaled_base is consistent with that when dealing with passing post-function-application representations.

May make more sense to use the regular scale operation and then "sanitize" it, with any negative radius replaced by a relevant change in angles.

Trying to ensure that radii remain positive, and taking particular care with UnitSphericalRepresentation.

mhvk · 2021-07-13T20:17:45Z

@adrn - would you be able to review this? Otherwise, perhaps fine to let @nstarman have a last look?

nstarman · 2021-07-13T21:55:47Z

    def __neg__(self):
-        self._raise_if_has_differentials('negation')
-        return self.__class__(self.lon + 180. * u.deg, -self.lat, copy=False)
+        if any(differential.base_representation is not self.__class__


I think a short comment would be good.

nstarman · 2021-07-13T22:09:22Z

+            in longitude should not be scaled if its spherical base is scaled
+            in radius.
        """
        scaled_attrs = [op(getattr(self, c), *args) for c in self.components]


Elsewhere you use generator comprehension, not list comprehension. I think the newer pattern is better.

Suggested change

scaled_attrs = [op(getattr(self, c), *args) for c in self.components]

scaled_attrs = (op(getattr(self, c), *args) for c in self.components)

nstarman

It looks good! but I would like @adrn to take a look as well. I'm still worried about the negative distances, but it doesn't seem to raise any errors, which seem comprehensive!

nstarman · 2021-07-13T22:21:02Z

+@pytest.mark.parametrize('op,args', [
+    (operator.neg, ()),
+    (operator.mul, (10.,))])
+def test_operation_unitspherical_with_rv_fails(op, args):


Is this the test for line 1667 in representation.py?

adrn

Sorry for the long delay here! I think this looks mostly great, but I had one question about the difference between UnitSphericalRepresentation and UnitSphericalDifferential. I see the following using your branch:

>>> coord.UnitSphericalRepresentation(40*u.deg, -11*u.deg) * 15*u.pc
<SphericalRepresentation (lon, lat, distance) in (deg, deg, pc)
    (40., -11., 15.)>
>>> coord.UnitSphericalDifferential(40*u.deg, -11*u.deg) * 15*u.km/u.s
<UnitSphericalDifferential (d_lon, d_lat) in deg km / s
    (600., -165.)>

I believe the TODO on L2051 points this out? But is there any reason to leave the asymmetry in the behavior here, or is it worth implementing the special-casing now?

mhvk · 2021-07-18T00:43:46Z

Sorry for the long delay here! I think this looks mostly great, but I had one question about the difference between UnitSphericalRepresentation and UnitSphericalDifferential. I see the following using your branch:
>>> coord.UnitSphericalRepresentation(40*u.deg, -11*u.deg) * 15*u.pc
<SphericalRepresentation (lon, lat, distance) in (deg, deg, pc)
    (40., -11., 15.)>
>>> coord.UnitSphericalDifferential(40*u.deg, -11*u.deg) * 15*u.km/u.s
<UnitSphericalDifferential (d_lon, d_lat) in deg km / s
    (600., -165.)>

The asymmetry here is because repr + derivative * time has to work as well, and for that case the multiplication of the derivative should just apply to all of its elements. This is somewhat annoying, and is the reason there is the scaled_base argument in _scale_operation -- so that if we do repr * scale, the unitsphericaldifferential is not multiplied by that scale (since the derivative of the angles is not affected by that scaling). Basically, differentials associated with representations do not work exactly the same as those which are not associated.

mhvk · 2021-07-19T18:50:18Z

@adrn - OK to merge this?

adrn

OK, that makes sense. Thanks @mhvk! Feel free to merge...

mhvk added coordinates Enhancement labels Apr 3, 2021

mhvk added this to the v4.3 milestone Apr 3, 2021

mhvk requested a review from adrn April 3, 2021 01:41

nstarman mentioned this pull request Apr 3, 2021

Add transform method to all representations #11444

Merged

pllim reviewed Apr 3, 2021

View reviewed changes

Comment thread docs/changes/coordinates/11469.bugfix.rst

mhvk force-pushed the representation-arithmetic-with-differentials branch from 73987cc to a6229e7 Compare April 3, 2021 17:03

mhvk marked this pull request as ready for review April 3, 2021 17:03

nstarman mentioned this pull request Apr 21, 2021

Transform radial rep when multiplication #11576

Merged

mhvk force-pushed the representation-arithmetic-with-differentials branch from a6229e7 to 50c1e27 Compare April 21, 2021 13:57

nstarman reviewed Apr 21, 2021

View reviewed changes

eteq modified the milestones: v4.3, v5.0 May 3, 2021

mhvk force-pushed the representation-arithmetic-with-differentials branch from 50c1e27 to 786b7d9 Compare May 4, 2021 00:03

nstarman reviewed May 4, 2021

View reviewed changes

mhvk force-pushed the representation-arithmetic-with-differentials branch from 786b7d9 to cd1bfdf Compare May 4, 2021 00:55

mhvk added 6 commits June 1, 2021 16:22

PEP8 changes

4f23bc6

Allow negative representation even if including differentials.

1052329

May make more sense to use the regular scale operation and then "sanitize" it, with any negative radius replaced by a relevant change in angles.

Allow scaling with differentials as well.

e8a6fa4

Trying to ensure that radii remain positive, and taking particular care with UnitSphericalRepresentation.

Changelog entry

35668a9

Adjust differential scale operation for base scaling

abf9fbd

Further small simplifications following review

b629e29

mhvk force-pushed the representation-arithmetic-with-differentials branch from cd1bfdf to b629e29 Compare June 1, 2021 20:24

mhvk mentioned this pull request Jul 13, 2021

Add registry of cosmology classes #11926

Merged

9 tasks done

nstarman reviewed Jul 13, 2021

View reviewed changes

nstarman approved these changes Jul 13, 2021

View reviewed changes

adrn reviewed Jul 17, 2021

View reviewed changes

adrn approved these changes Jul 19, 2021

View reviewed changes

mhvk merged commit 7063a8f into astropy:main Jul 19, 2021

mhvk deleted the representation-arithmetic-with-differentials branch July 19, 2021 20:29

		# super() checks that the class is identical so can this even happen?
		# (same class, different differentials ?)

	scaled_attrs = [op(getattr(self, c), *args) for c in self.components]
	scaled_attrs = (op(getattr(self, c), *args) for c in self.components)

Conversation

mhvk commented Apr 3, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions Bot commented Apr 3, 2021

Uh oh!

Uh oh!

mhvk commented Apr 21, 2021

Uh oh!

nstarman Apr 21, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nstarman Apr 21, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nstarman May 4, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk May 3, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

nstarman commented Apr 21, 2021

Uh oh!

eteq commented May 3, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

nstarman commented Jun 1, 2021

Uh oh!

mhvk commented Jun 1, 2021

Uh oh!

nstarman commented Jun 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nstarman commented Jun 1, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented Jul 13, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nstarman left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk commented Apr 3, 2021 •

edited

Loading

nstarman Apr 21, 2021 •

edited

Loading

nstarman Apr 21, 2021 •

edited

Loading

nstarman May 4, 2021 •

edited

Loading

mhvk May 3, 2021 •

edited

Loading

nstarman commented Jun 1, 2021 •

edited

Loading

nstarman commented Jun 1, 2021 •

edited

Loading