Sometimes angle units are dimensionless - new system, more general equivalency? - #1161

I think this might be worth bringing up on the mailing list.

I wonder if rather than worrying about different 'systems' or 'equivalencies' it should be possible to have units of "rad" that are dimensionless and can be freely dropped, but are otherwise included in a quantity until/unless they're dropped.

Or maybe u.rad could be callable. For example u.rad(dimensionless=False) would return an irreducible version of radians. I think by default it should be dimensionless--I don't usually care to have "rad" showing up in my units. But I agree there are times it can be helpful, especially for bookkeeping.

It doesn't appear this generated much interest on the list, so I'll just comment here.

I don't really see why it ever makes sense to drop it implicitly - even though an angle is formally dimensionless mathematically, I'm not sure there's any context where it doesn't help to attach the "angle" label to it. Any sort of implicit dropping is bound to get people into trouble, because what might be the "obvious" place to do it for one person will probably seem unnatural to another with a different background.

I don't entirely follow the example you gave @mhvk - in particular, I'm not sure what P is supposed to be (a Quantity with time units, or a raw number?), but it seems to me you could either do Erot = 0.5 * I.to(u.erg) * Omega**2 or Erot = 0.5*I*Omega.value**2 ; Erot.to(u.erg) and get out what you want.

The one case I do see value here is that right now we want np.sin(someanglequantity) to give out a dimensionless quantity, but that will hopefully be addressed in #929 or some related follow-on.

Oh, and regarding the comment on-list about having rad be derived from dimensionless: the problem there is then it doesn't have a physical_type of angle, which I would think we want.

@eteq - Yes, it doesn't seem to have been picked up. I also mentioned the issue in #1139 - let me ping @astrofrog, @mdboom, and @iguananaut here, so perhaps they can comment.

On the nature of the angle: SI defines radian as a unit derived from dimensionless, so in that sense there is no ambiguity. That said, I do see the advantage of having an angle physical type in catching silly mistakes.

About my example: sorry I was not clearer; indeed P is a period in units of time, Omega the associated angular frequency, and I the moment of inertia. For a given I, Omega, the standard expression for the rotational energy is the one I list, Erot = 0.5*I*Omega**2. However, if rad is a unit of angle, then Erot cannot be converted to erg:

import astropy.units as u
import numpy as np
P = 2.*u.s
Omega = 360. * u.deg / P
phase = Omega * 1.* u.hr
I = 1.4 * u.M_sun * (10*u.km)**2
Erot = 0.5*I*Omega**2
Omega, phase, np.sin(phase), Erot

gives

(<Quantity 180.0 deg / s>,
 <Quantity 180.0 deg h / s>,
 <Quantity -6.11403104275e-13 >,
 <Quantity 2268000.0 deg2 km2 solMass / s2>)

which is all correct. But

Erot.to(u.erg)
ERROR: UnitsException: 'deg2 km2 solMass / s2' and 'erg' (energy) are not convertible [astropy.units.core]

As it is now, it is not easy to do this elegantly, except by writing explicitly (Omega.to(u.rad/u.s).value * 1./u.s).
My question is whether we can implement a rad=unitless system, so I could do something like Erot.to(rad=1), or Erot.to_system(rad=1). Or perhaps just treat angles as dimensionless, following SI, and live with some risk of people taking np.sin of a number that is dimensionless but makes no sense as an angle.

I'm not really comfortable with having two different systems, so if we decided to make rad dimensionless, it would need to always be dimensionless. One option is also to simply not use it when one wants to use it in a dimensionless way. I'd be interested to hear what other people think.

Ah, I understand the example now @mhvk (the M_sun was throwing me off...). I now better understand the issue. So with that in mind, it's not clear to me that there's really any perfect solution to this problem, because angles do have a physical meaning distinct from some other dimensionless numbers... But this is also a clear example of how sometimes they shouldn't be. So no matter what we do it's going to lead to ambiguity in some cases.

So that leads to the suggestion of having some quick but explicit way of converting between radians and dimensionless units. I'm not a big fan of Erot.to_system(rad=1), though, because that's a major change to the signature of to (which currently must accept a Unit or something that can be made into a Unit. And while I see what you're doing with that, my experience is that using syntax "tricks" like that tend to lead to substantial confusion among novice users (or even experienced ones who don't read code too carefully).

So I tend towards agreeing with @astrofrog second suggestion that we should just require that users force it to be dimensionless when they want to use it this way. In your example's case, that would probably most-mean doing Omega = (360 * u.deg).to(u.radians).value / P or something like that (or Omega = 2 * pi / P, but I admit that's cheating if you're starting with an angle from somewhere else).

But I agree with @astrofrog it would be useful to hear what others think of that approach.

Oh and one other important detail that just occurred to me: we really do need angle to be a distinct physical type, because converting the Angle class to being a subtype of Quantity (important for a variety of things in the coordinates package), depends on angles being distinct from other dimensionless quantities.

@eteq, @astrofrog - I certainly see the arguments in favour of having angle as a real dimension, but it is an explicit decision not to be consistent with the current definition of SI [1]. Do we really want that? Though one can certainly argue we follow good precedent here, as the IAU still lists angles as "supplementary units" [2], which is how SI used to treat them.

In any case, it seems we agree there should be an option. My sense is that it should be possible to be consistent with "strict" SI rather than going via .value (which is a hack that I feel should be reserved for cases where "it is unphysical, but I know what I'm doing"). The main question would be one of implementation. Three options would seem to be

1. a si_strict system, in which rad and sr really are derived units (or, equivalently, make the current .si follow SI as currently defined, but introduce an .si_with_angles or .si_old or .iau or so that keeps rad as a unit in its own right). This avoids polluting .to but adds another system beyond the existing .si and .cgs.
2. an option to .to similar in nature to equivalencies, but broader, i.e., "if possible, substitute as many rad for 1 and vice versa to get the conversion to work"); essentially, it would remove all rad in the current and target units, convert as appropriate, and then add target rads back.
3. An alternative to .to that does the above.

Of the three, I prefer the first option (with .si being current SI; old SI having a different name).

[1] http://physics.nist.gov/cuu/Units/units.html
[2] http://www.iau.org/science/publications/proceedings_rules/units/

I was actually going to suggest solving this with a built-in equivalency (a la suggestion (2)) above. I think that's probably the cleanest solution and will make everything work. The problem is that the equivalencies don't currently do any factoring, so it won't find the match in the compound unit you gave in your example. I had meant to add that once the compose functionality was fully working, I just haven't got around to it yet. So let's keep this issue open as a reminder to get that done.

@mdboom - OK, sounds good. Out of curiosity, if one were to implement a strict SI system, how would one actually define a unit derived from the dimensionless one? I tried something probably insensible (certainly didn't work...):

def_unit(['rad', 'radian'], 1 * Unit(1), register=True, prefixes=True,
         doc="radian: angular measurement of the ratio between the length "
         "on an arc and its radius")

I'm slightly puzzled, since I can do percent = Unit(0.01) (although this of course doesn't keep the name)

@mhvk: In our system, all units are (logically) derived from dimensionless by default. There isn't really a special token for dimensionless, it's just indicated by the lack of units. Whether that needs to be re-addressed isn't clear -- I'll look into it when I implement a solution to this.

@mdboom - still leaves me a bit confused: I can set 1 * u.deg and it will keep it as deg even after combining with other units, until I explicitly ask to convert to rad. But I guess here rad is the irreducible unit, and perhaps there is no equivalent for dimensionless?
It would be handy more generally to have (the ability to define) units that are dimensionless; apart from % alluded to above, one could have the mma (= milli-modulation amplitudes) typically used for white-dwarf pulsations.
A class more similar to rad is of course the infamous magnitude, which really is dimensionless, but is also a unit that is often combined with others (transient decay rates are often measured in mag/day. (But let me stop before I drive you to despair... ;-)

Right. In our system we essentially equate dimensionless with having no unit. And it's the result of dividing any unit by itself (since the units cancel out and we are left with no units).

In [2]: u.m / u.m
Out[2]: Unit(dimensionless)

So dimensionless is not actually a unit from which one can derive, but a special case meaning "has no units".

I believe what you're suggesting is having a dimensionless "unit" from which some other units derive. But that would break the very useful property above, would it not?

To elaborate further, since dimensionless is the special case meaning "having no units", I'm not sure why it's helpful to have units derived from it, other than perhaps having a little syntactic sugar that's completely ignored by the system.

Barring a small buglet that can be fixed thusly (see diff below), it is possible to define units in terms of dimensionless, but you don't really gain anything, since any dimensionless unit would than be convertible to any other. The problem with that is illustrated by then deriving rad from it. It then means you can't distinguish between then angle, solid angle (since 1 ** 2 == 1), angular momentum from momentum etc. It seems to just punch a large hole through what this system is trying to do by throwing out information. What exactly would you be trying to achieve by defining named dimensionless units (other than making the code more self-documenting)? Is there anything the unit system can and should enforce with these things?

diff --git a/astropy/units/utils.py b/astropy/units/utils.py
index b620ab3..e0a76b5 100644
--- a/astropy/units/utils.py
+++ b/astropy/units/utils.py
@@ -58,7 +58,11 @@ def generate_unit_summary(namespace):
             continue

         if isinstance(val, core.PrefixUnit):
-            has_prefixes.add(val.decompose().bases[0].name)
+            decomposed = val.decompose()
+            if len(decomposed.bases):
+                has_prefixes.add(val.decompose().bases[0].name)
+            else:
+                has_prefixes.add('dimensionless')
         else:
             units.append(val)

Yes, it would just be syntactic sugar, but with an "SI approved" stamp! Note that the NIST site defines
plane angle rad = m * m**(-1) = 1, so your example is appropriate ;-).

Kidding aside, your change makes sense for % and mma at least: those are just similar things to having the ability to convert AU to m (if with more boring conversion factors). So, in that sense, I would argue to go ahead and implement it.

The question whether rad should also be derived from it, is a separate one. I think it all boils down to well-known problems with angles, which apparently SI couldn't really make their mind up about either (since they changed from "supplementary" to "derived" units) and hence we probably cannot solve it either. But I think with your small change it becomes easy to offer people the option to choose what is most convenient (whether with the now very easy different "system" option, or still via an equivalency).

👍 to @mdboom's solution of using a built-in equivalency. Now that I think about it, that's exactly the whole point of equivalencies! (situations where you want to convert units in special/context-dependent ways.)

I'm assigning myself to this. If I can find a reasonable solution to this in time for 0.3.0, great, otherwise I don't consider it highly urgent.

Agreed it is not hugely urgent, though in the meantime I did encounter another fairly common situation where it would help if radian was equivalent to dimensionless:

exp(1j*angle)
exp(1j*angularfrequency*time)

(and likely other examples with complex numbers).

@mdboom - I made a trial implementation. It is somewhat ugly, since really this is not quite the same as a normal equivalency, in that I'd want all powers of radian to be converted, independent of whether it is part of a more complex unit. For equivalencies like u.spectral(), that is not the case -- one really does not want kg*m/s^2 to become kg/s^3 is I use u.spectral()!

Ideally, I'd like this to be an option to the user, i.e., they can set somewhere at the top of the script whether they want this to be the case. I now wrote this as an attribute to _current_unit_registry, but you can probably think of much smarter ways

Since it is a bit of a weird equivalency, I also tried to write a PR where this is by option only; see https://github.com/mhvk/astropy/tree/units-radian-is-unity-decompose

Anyway, with this, I can now do:

import astropy.units as u
Omega = u.cycle/(1.*u.s)
I = 1.4 * u.M_sun * (10*u.km)**2
Erot = 0.5*I*Omega**2
Erot.to(u.erg, equivalencies=u.radian_is_dimensionless())
# or  Erot.to(u.erg, equivalencies=((u.rad,None),))

which yields

<Quantity 5.496856431977914e+46 erg>

Yeah -- it's a bit hackish, but I see why it has to be the way it is. I think I'm okay with this.

This will need some tests before merging.

@mdboom - OK, I can write some tests. What I'd like your opinion on, however, is how to set this script/context wide. If that is not immediately obvious, perhaps I'll cut the current lines that do:

        if getattr(_current_unit_registry, 'radian_is_dimensionless', False):
            normalized += [(Unit('rad'), None, lambda x: x, lambda x: x)]

and we can do a separate PR on that.

Well, there's at least three choices here:

1. We have global equivalencies that are stored somewhere. These would always add themselves to whatever equivalencies were explicitly passed in. If we go this way, I would just put these in the current registry in its own list of equivalencies, not special casing it like you've done here.

2. We have equivalencies as a member on Quantities. I'm not a fan of that, but I know we do that now. But that doesn't help with "bare" units.

3. Neither, and equivalencies are always explicitly passed. I actually prefer this option because it avoids global state, object state or other magic that could make things confusingly behave differently based on context.

I feel somewhat strongly about being able to use "true SI" in Quantities/Units without having to add the equivalency everytime the quantity is initialised or converted. Do you think your context manager for imperial units be repurposed for this? I.e., something like

from astropy.units import radian_is_dimensionless:
radian_is_dimensionless.enable()

or

with add_enabled_units_context(radian_is_dimensionless):
    ...

Though perhaps I should not think of repurposing. Would it be easy to add something like:

with add_default_equivalency(u.radian_is_dimensionless()):
    ....

(this may also be an alternative to storing equivalencies in quantities...)

Your last suggestion is pretty close to what I meant by my option (1) above, but I'm still not crazy about introducing more global state.

@astrofrog - it appears the travis failure to build sphinx is due to modeling? At least, assuming it is due to the RunTimeWarnings from that package.

@astrofrog, I raised an issue about what the warning I saw, though I'm no longer sure it is related, as I also saw it in #1593. The latter also has other problems, so I'll at least fix those.

Yay! it is in good standing again. I think this is ready to be merged; @mdboom, @astrofrog?

This seems fine to me, but I think that dimensionless_angles would be much more readable than angles_dimensionless - what do you think @mhvk?

@astrofrog - I actually had that initially, as I agree it reads better. The main reason I changed is that I use dimensionless_unscaled a lot, and would love to keep it easy to TAB-complete that...

An alternative would be to go back to radian_is_dimensionless, but I think it is better to have the name include the physical item.

I can see your point about interactive usage, but really, I don't think that should determine the wording, otherwise as time goes on we'll add convoluted terminology just so that it tab completes well ;) I personally think dimensionless_angles or dimensionless_radian would be the clearest.

@astrofrog - OK, it is back to dimensionless_angles...

This looks ready to merge.

However (and apologies for this), we seem to be working on top of each other a bit today 😉. In light of #1662, do you think you could rebase this on #1662 and make the analogous changes so there's no set_enabled_equivalencies/set_enabled_equivalencies_context dichotomy? We can always fix this up after the fact of course, but I'd rather not miss it -- I think #1662 makes things much simpler/cleaner/better.

@mdboom - the only reference to the context manager in this PR was a test marked xfail; this should not have been here in the first place! Anyway, I now removed it, which should make this PR at least independent of #1662.

Finally, it seems all is set here.

@mhvk - it looks ok and Travis passes. It sounds like #1662 should be merged first? (@mdboom?)

Yeah, I'll merge #1662 first, once it gets through Travis, then this. #1662 may force a rebase of this -- I'm happy to manually merge if necessary. We'll see.

@mhvk - this needs a rebase now that #1662 is in - does it look like it's going to be complicated?

I think I can take care of the rebase... hold on...

Merged (manually)

@mdboom - thanks! Let me know if #1593 is a problem.