Consider whether bounds are needed as part of the Angle class. - #1310

As discussed in #1006, there are problems that arise when the bounds of an Angle are explicitly enforced, particularly with regard to addition and subtraction. The current behavior is that two Angles can be added or subtracted if they have the same bounds, and the output will have the same bounds.

An example problem is the Haversine distance, which performs a subtraction of dec2 - dec1. This can end up having any value between -180 to +180 (versus the latitude bounds of (-90, 90)). Therefore the Angle class cannot reasonably enforce bounds by raising an exception when an algebraic result is outside the bounds of the input values. To do so would break a lot of code that should just work, so it seems that raising exceptions for out-of-bounds values is undesirable.

For truly circular angles the bounds specify where to wrap. This can be useful, but can lead to problems for expressions like sin((a + b) / 2). For a=359 and b=3 this will return (with wrapping) the value sin(1 degree) which is different from the expected sin(181). So doing wrapping has the potential to produce numerically incorrect results for trigonometric calculations with no warning. This would imply that wrapping is dangerous and should be disabled unless arithmetic and ufuncs are disabled.

Without wrapping or exceptions, then the bounds are not really doing much as far as I can tell.

cc: @adrn @demitri @mdboom @eteq

I would argue that for subtraction, the result should no longer have any bounds for the reason you mention, so if we want to keep bounds, we certainly have to improve the rules on when they are kept and when they are removed.

Subtraction is just another name for addition (a + -b), but if we remove bounds from the result of addition or subtraction (and negation, multiplication, and division) that will address many of the issues.

I think this defeats the purpose of bounds. They should not be silently dropped.

I think this defeats the purpose of bounds. They should not be silently dropped.

@demitri - Do you have any suggestion for handling the issues that I mentioned?

@taldcroft Right now, no. I'd have to sit down and analyze the problem rather than answering off the cuff. I'd say disable the operations under question for the moment.

My main guiding thoughts:

• Bounds should not be dropped on operations. One purpose of them is to be a "safety net", and as such should not be removed. I'd rather the code fail or bounds not implemented (which is a shame, as I think there's read value there).
• Note that there is a difference between Angle and Coordinate. If the latter needs some special treatment (e.g. the Haversine distance is an issue with lat/lon, not bare angles), then that behaviour belongs in the subclass.
• Note that things like the means of angles are not linear operations (see above).

I think this would be worth presenting as a discussion on the astopy list.

Disabling basic arithmetic and ufunc operations would certainly spoil the joy of having merged #1006, which was a lot of work! Here is what we currently get, which I find unacceptable (and will open an issue soon):

from astropy.coordinates import ICRSCoordinates
import astropy.units as u
c1 = ICRSCoordinates(0, 89, units=(u.deg, u.deg))
c2 = ICRSCoordinates(0, -89, unit=(u.deg, u.deg))
from astropy.coordinates.angle_utilities import haversine_sphere_dist
haversine_sphere_dist(c1.ra, c1.dec, c2.ra, c2.dec)

ERROR: BoundsError: The angle(s) -178.0 falls outside of the specified bounds (-90.0, 90.0) (in deg) [astropy.coordinates.angles]
---------------------------------------------------------------------------
BoundsError                               Traceback (most recent call last)
<ipython-input-8-2452654b2f4b> in <module>()
----> 1 haversine_sphere_dist(c1.ra, c1.dec, c2.ra, c2.dec)

/data/baffin/tom/git/astropy/astropy/coordinates/angle_utilities.py in haversine_sphere_dist(lon1, lat1, lon2, lat2)
    607     from math import asin, sin, cos
    608 
--> 609     sdlat = sin((lat2 - lat1) / 2)
    610     sdlon = sin((lon2 - lon1) / 2)
    611     coslats = cos(lat1) * cos(lat2)

/data/baffin/tom/git/astropy/astropy/coordinates/angles.py in __sub__(self, other)
    280     def __sub__(self, other):
    281         self._bounds_match_check(other, 'subtract')
--> 282         return super(Angle, self).__sub__(other)
    283 
    284     def __mul__(self, other):

/data/baffin/tom/git/astropy/astropy/coordinates/angles.py in __array_wrap__(self, obj, context)
    260 
    261         if isinstance(obj, Angle):
--> 262             return Angle(obj.value, obj.unit, bounds=obj.bounds)
    263 
    264         return obj

/data/baffin/tom/git/astropy/astropy/coordinates/angles.py in __new__(cls, angle, unit, bounds, dtype, equivalencies)
    133             raise u.UnitsException("No unit was specified")
    134 
--> 135         bounded_angle, bounds = cls._bounds_check(angle, bounds, unit)
    136 
    137         self = super(Angle, cls).__new__(

/data/baffin/tom/git/astropy/astropy/coordinates/angles.py in _bounds_check(angle, bounds, unit)
    224                 angle = np.where(too_small, angle + TWOPI, angle)
    225                 if np.any(angle > upper_angle):
--> 226                     raise_error(original_angle, lower_angle, upper_angle, unit)
    227 
    228                 too_small = angle < lower_angle

/data/baffin/tom/git/astropy/astropy/coordinates/angles.py in raise_error(original_angle, lower_angle, upper_angle, unit)
    176                     u.radian.to(unit, lower_angle),
    177                     u.radian.to(unit, upper_angle),
--> 178                     unit))
    179 
    180         if bounds is None:

BoundsError: The angle(s) -178.0 falls outside of the specified bounds (-90.0, 90.0) (in deg)

@taldcroft This is what I mean about Angle not being Coordinate. Maybe Angle bounds should be dropped when they are used as components of coordinates, and the Coordinate enforces its own bounds (e.g. dec in [-90,90]).

@taldcroft, @demitri - Just read the discussion in #1006 and here. My feeling is that the bounds in Angle do not work in part because one cannot easily capture both a latitude and a longitude in one class. Indeed, arithmetic on one is not possible without knowing the other (e.g., the average of 30 and 30 is 30 for a longitude, but usually not for a latitude; and the average of longitudes 30 and 31 does depend on the corresponding latitudes).

To me, this suggests bounds checking should not be part of Angle. But it could still happen in two different places:

(1) In RA, Dec (which, however, are just specific examples of more general Longitude and Latitude classes used generally by coordinates). For these, bounds are checked and only a few selected operations are enabled, a bit like Time. E.g., sum/division/multiplication/mean/etc. are not possible, subtraction returns an Angle (closest to zero for longitude), addition/subtraction of an angle is possible (bound check for latitude, wrap for longitude), comparison is allowed (obvious for latitude; using "shortest distance" for longitude, i.e., 23h < 1h is True, because 23 is closer to 1 on the low side).

(2) Bounds checking becomes part of the coordinate class (if so, RA and Dec should probably be dropped).

As I think @astrofrog mentioned, in practice we may be beating a dead horse, as for use in coordinates one does not really need to use angles directly at all. But I do see some use outside of coordinates, mostly as a neat way to represent/parse angle quantities; for those, however, I think bound checking just confuses things.

@mhvk - I agree that bounds (as a property that either enforces limits or does auto-wrapping) should be dropped from Angle.

If I understand correctly, your option (1) would mean that for a Coordinate object c then c.ra and c.dec would be objects that don't support addition, multiplication or division. I think that would be pretty limiting, and in particular would break a lot of standard spherical trig computations, including the astropy spherical distance functions. That seems like a show-stopper to me.

That leaves option (2), which seems good to me. I personally don't see a real need for RA and Dec classes, this just seems to add unnecessary complication. If you have RA and Dec then why not also Alt and Az and L and B classes?

The Angle class does provide important utility in terms of formatting in the arcane systems that astronomers use. This is a very good thing for c.ra, c.dec, or c.b. After dropping bounds there is one thing missing which is a way to specify the output representation in the range (-180, 180) vs. (0, 360) or in general whatever 360 slice you want.

Just quickly brainstorming I was thinking I was thinking about a wrap_angle attribute that specifies where in the circle to wrap when representing the angle as a string. There might also be a method that can return a new Angle with the internal values wrapped accordingly, but importantly the action of wrapping would only be when the user explicitly requests it, or when formatting as a string. Note that there is no explicit support for a latitude range of (-90, 90), and I don't think this is ever needed, but counter arguments here are welcome!

Without a bounds attribute then the mean method really needs to be just the plain old linear mean instead of circular quantity mean (http://en.wikipedia.org/wiki/Mean_of_circular_quantities) because there is no way to uniquely assign the output phase. I.e. what is the circular mean of [721, 1081]? Should it be 1, 721, 1081 or none of the above? I think it is less surprising if the mean method is just the linear mean. However we could introduce a new method like circular_mean which is very obviously the mean as defined for a circular quantity,with a well-documented output range of [0, 360].

@taldcroft - Indeed, the coordinates class would be composed of a longitude and latitude that have restricted ability. I don't think this is necessarily a problem: for transformations, all one needs is to be able to subtract angles (which would still work, just returns an Angle rather than a longitude/latitude) and take trigonometric functions. But in the end this seems more a detail of how one implements coordinates; the only real advantage I can see for having the separate subclasses is that the longitude one would also be useful for a polar coordinate system.

As for bounds, it would seem that one could use the existing mechanism in the longitude/latitude classes (or the equivalent inside coordinates).

As for np.mean, I think it is a risk operation for Angle, but see no reason to disallow it. For longitude/latitude, it is not uniquely defined on the individual angle, so should not be possible (though there is no reasons one could not average coordinates).

It seems to me that no one is really treating the Angle object as a pure geometric object representing an angle. Nearly all discussion about it pulls in coordinates, which it should not: Angle is a base class. If the line is blurred here, then I would just propose making Angle private since (to my mind) the implementation is not what I'd think it would be. Then it can be made exclusively to support coordinate objects with wild abandon.

@taldcroft The idea of a wrap_angle overlaps the concept of bounds far too much - I think in the end it's a slightly different name. The idea of having RA and Dec classes was to provide convenience objects since the vast majority of coordinate representations are ra and dec. They'd have the expected bounds checking, etc. There's no need to create the Alt, Az, etc. since they are not nearly as common, and L and B as class names looks silly. It's not like all coordinate systems are implemented in roughly equal proportions: the overwhelming majority are RA/Dec, so these classes are just convenient.

@demitri - I entirely agree that Angle is a pure geometric object, which is part of the reason I'm arguing to drop bounds. For an object which should support generic trig computations (which tend to involve add, subtract, divide, multiply), having a bounds attribute breaks everything. It just doesn't work. Ideas welcome of course.

It happens that our main use case for Angle is to represent each of the two components in a coordinate, that's why alot of the discussion has been based on that context. We want to be sure that Angle works perfectly in that case.

wrap_angle is meant to provide limited wrapping functionality within an Angle object that has no bounds attribute, so there would be no overlap.

Just looking at the actual code for RA and Dec, they basically do nothing except put in bounds. (OK, RA has a couple of trivial convenience methods that could easily be standalone utility functions.) I would propose a simplification by just eliminating RA and Dec as subclasses since they aren't doing anything. This will let us cut a bunch of boilerplate code and docs [1]. If RA and Dec are purely providing a trivial syntactic sugar, then they could be dropped.

[1] Kenneth Rietz at PyCon: "My favorite pull requests are ones that delete code".

As a proof of concept I just spend 15 minutes hacking bounds out of the Angle class, in the process removing 225 lines of code and making a number of simplifications. The only tests I removed were ones that were specifically checking bounds functionality. With one exception (test of hour_angle method), the all tests passed on the first try. I think this is a pretty good indication that bounds are not needed in Angle and removing it just makes everything simpler.

As an immediate benefit, all our own distance functions like haversine_sphere_dist actually work with Angle! [1]

Of course there is a caveat that coordinates test coverage might not be complete with respect to wrapping issues, since this was known to be handled previously within Angle.

[1] OK, after changing math trig functions to numpy functions, as will be done when coords is array-ified.

@taldcroft - just an idea - how about reinstating the bounds as an option in the initializer (maybe under a different name) to allow the angle to be automatically wrapped? (would prevent the user from having to do the arithmetic themselves).

I think this discussion should be taken to the astropy-dev list. @taldcroft, since you originally suggested this change, do you want to do that?

Bounds are useful for that reason alone - to not have to write your own code to wrap results. No need to change the name - that's the purpose. That and providing a way to restrict values for error checking.

What about:

• Initializer takes a bounds argument for wrapping, but then discards the bounds after wrapping the value in the initializer
• Angle has a wrap method that can be used to wrap the value to a set of bounds, raising an exception if the value is outside.

Then this allows arbitrary arithmetic between Angle objects, but provides a facility for wrapping and raising exceptions if needed after an operation.

@astrofrog How could/would this use case be handled? I define declination:

dec = Angle("80d", ???)
dec + Angle("50d")
<Exception: This angle cannot exceed 90°.>

Obviously I don't want the result to be 30°.

@demitri - well in my scenario in Angle you could do:

>>> dec = Angle("80d", bounds=(-90, 90))
>>> dec2 = dec + Angle("50d")
>>> dec2.wrap(-90., 90)
<Exception: This angle cannot exceed 90°.>

Of course, Angle is not a declination so this would not be available by default. The Declination class on the other hand could support raising an exception using this mechanism behind the scenes.

@astrofrog That usage of wrap defeats the purpose of bounds checking if I have to manually perform the check every time the angle might exceed my stated constraints. That's the whole point of building it into the class.

No, Angle is not a declination, but I should be able to define an Angle object with arbitrary bounds. And there was discussion of removing the Dec and RA classes, something I disagree with for exactly these reasons. Why should everyone have to implement the very commonly used declination object?

@demitri - in my scenario the purpose of wrap is to be available for sub-classes such as Declination, which exists in astropy.coordinates. So you are incorrect to say that everyone would have to implement it - we can provide both Angle and Declination. Angle would not do bound-checking but would offer methods that can do it, and Declination would make use of it.

@astrofrog I was making an argument for not dropping Dec as a class. Does this mean if I want the same behaviour with different bounds (e.g. [-45°, 45°]) that I would have to make a subclass? This is what seems unnecessary and should be supported by Angle natively.

Just out of interest, can you give a real life example where the bounds are not -90 to 90, -180 to 180, or 0 to 360? It would help me understand why we need arbitrary flexibility in this respect.

I don't have a specific example, no. It just seems to me that if we support one, it should be a similar effort to support any. Well, wait. Let's say someone was working on coordinates for a telescope pointing. What if, due to geometry, mechanics, trees, whatever, the telescope could only point to targets with a declination range of [26.7, 82.1]? Not unreasonable, and it would be really nice to have code that made sure that the angle I give it falls in that range without me having to check all the time. I could very easily forget to check.

Another use case is using the Angle class for rotations. I could define an angle of 400° to mean than it's one full rotation plus 40°, but wouldn't want that value wrapped to be 40°.

For the latter example, that's just the standard behavior in Angle though, which we aren't suggesting changing.

What about:
Initializer takes a bounds argument for wrapping, but then discards the bounds after wrapping the value in the initializer
Angle has a wrap method that can be used to wrap the value to a set of bounds, raising an exception if the value is outside.

Then this allows arbitrary arithmetic between Angle objects, but provides a facility for wrapping and raising exceptions if needed after an operation.

@astrofrog - this is pretty reasonable, but I think it can be even simpler and meet all typical use cases. What about just the following two methods and nothing cluttering up the initializer. These two methods explicitly decouple wrapping and bounds checking:

• wrap_at(angle): wraps within 360 degrees centered on the specified angle. So to get -180 to 180 you use wrap_at('180d'). To get -270 to 90 you wrap_at('90d'). To ensure 0 to 360 you wrap_at('360d') or wrap_at(360 * u.deg).
• check_bounds(angle1, angle2): check that all angles are within angle1 to angle2 and raise an exception if not.

Some typical code would be:

ra180 = ra.wrap_at(180 * u.deg)  # convert from 0-360 to -180 to 180
dec.bounds_check(-90 * u.deg, 90 * u.deg)  # validate our dec values

This makes for explicit code that is readable and allows intuitive math and spherical trig with Angle objects.

I think I'm ok with either of @astrofrog or @taldcroft's suggestions here, but one thing is still left in my mind: what about the string representation? There are a lot of cases where you want the string representation to also be wrapped to a certain range (e.g. 0 to 24 hr for RA, -90 to 90 deg for dec). Right now that's implicitly taken care of by the bounds, but if that's removed, then subclasses that want this would have to re-implement string formatting code, right?

Looks good to me! (apart from the docs needing updating)

OK, I believe I have addressed all comments and done all the requisite doc updates. From my perspective this is good for merge.

The items above are all pretty straightforward (except the array-output question, which is probably for another PR). So aside from that, I think this is ready to merge. Thanks for all the work, @taldcroft !