Alternative astropy.constants API - #572

I still kind of feel iffy about the inconsistency of this... I would still prefer to be able to from astropy.constants import e. Since this already adds a special EMConstant class (which seems reasonable), why not go ahead and add EMConstant.esu and EMConstant.gauss properties--update the error messages accordingly? The e instance by itself isn't usable in combination with other quantities, but e.esu and e.gauss would be. I'm not exactly sure how things work when combining constants with other quantities in other units, but I think something could also be kluged together so that e doesn't combine with SI-compatible units without using e.si explicitly too.

I'll see if I can some up with something illustrating what I'm thinking, if it will even work.

By the way, I'm getting two test failures in the cosmology package with this PR:

================================================= FAILURES =================================================
_______________________________________________ test_ogamma ________________________________________________

    @pytest.mark.skipif('not HAS_SCIPY')
    def test_ogamma():
        """Tests the effects of changing the temperature of the CMB"""

        #Tested against Ned Wright's advanced cosmology calculator,
        # Sep 7 2012.  The accuracy of our comparision is limited by
        # how many digits it outputs, which limits our test to about
        # 0.2% accuracy.  The NWACC does not allow one
        # to change the number of nuetrino species, fixing that at 3.
        # Also, inspection of the NWACC code shows it uses inaccurate
        # constants at the 0.2% level (specifically, a_B),
        # so we shouldn't expect to match it that well. The integral is
        # also done rather crudely.  Therefore, we should not expect
        # the NWACC to be accurate to better than about 0.5%, which is
        # unfortunate, but reflects a problem with it rather than this code.
        # More accurate tests below using Mathematica
        z = np.array([1.0, 10.0, 500.0, 1000.0])
        cosmo = core.FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=0, Neff=3)
        assert np.allclose(cosmo.angular_diameter_distance(z),
                           [1651.9, 858.2, 26.855, 13.642], rtol=5e-4)
        cosmo = core.FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725, Neff=3)
>       assert np.allclose(cosmo.angular_diameter_distance(z),
                           [1651.8, 857.9, 26.767, 13.582], rtol=5e-4)
E       assert <function allclose at 0xa064320>(array([ 1651.91440293,   858.20451145,    26.85856246,    13.64688241]), [1651.8, 857.9, 26.767, 13.582], rtol=0.0005)
E        +  where <function allclose at 0xa064320> = np.allclose
E        +  and   array([ 1651.91440293,   858.20451145,    26.85856246,    13.64688241]) = <bound method FlatLambdaCDM.angular_diameter_distance of FlatLambdaCDM(H0=70, Om0=0.3, Ode0=0.7)>(array([    1.,    10.,   500.,  1000.]))
E        +    where <bound method FlatLambdaCDM.angular_diameter_distance of FlatLambdaCDM(H0=70, Om0=0.3, Ode0=0.7)> = FlatLambdaCDM(H0=70, Om0=0.3, Ode0=0.7).angular_diameter_distance

astropy/cosmology/tests/test_cosmology.py:113: AssertionError
___________________________________ test_angular_diameter_distance_z1z2 ____________________________________

    @pytest.mark.skipif('not HAS_SCIPY')
    def test_angular_diameter_distance_z1z2():
        tcos = core.FlatLambdaCDM(70.4,0.272,Tcmb0=0.0)
>       assert (tcos.angular_diameter_distance_z1z2(1, 2) ==
                646.22968662822018)
E       assert 646.22968662822007 == 646.2296866282202
E        +  where 646.22968662822007 = <bound method FlatLambdaCDM.angular_diameter_distance_z1z2 of FlatLambdaCDM(H0=70.4, Om0=0.272, Ode0=0.728)>(1, 2)
E        +    where <bound method FlatLambdaCDM.angular_diameter_distance_z1z2 of FlatLambdaCDM(H0=70.4, Om0=0.272, Ode0=0.728)> = FlatLambdaCDM(H0=70.4, Om0=0.272, Ode0=0.728).angular_diameter_distance_z1z2

astropy/cosmology/tests/test_cosmology.py:390: AssertionError

I'm on scipy 0.10.1. The second one looks like just a slight fp error, and should be using assert_close. The other I'm not sure about I need to take a closer look.

+1 to this w/ @iguananaut's suggestion of having these work as e.si and so on.

While this allows usage like e.esu.si, I think that's actually good, because that actually says "use the ESU convention for what e means, but then cast it into the cgs unit of charge. This would be equivalent to e.esu.to(u.coulomb) which works in both cases.

I've included a commit that @iguananaut wrote which implements his suggestions. I agree with his comment in the commit message that e on its own should probably not be usable as a quantity in arithmetic operations without specifying e.si.

@iguananaut - do you want to try and implement the fix to make it so e cannot be used without specifying the system? I started to try and overload __add__ etc for EMConstant but that doesn't work because then e.gauss can also not be used for arithmetic. Maybe we need a third class, EMConstantGroup, that does not inherit from Quantity and can show the available systems when doing __repr__ or _str__?

@astrofrog Thanks for pulling my suggestions--I'm glad it seems to work out. As for overriding the operators, I think it would make sense for EMConstant instances to have some built in knowledge of what system they belong to. They should have a .system attribute and a system=None argument in their initializer. If .system is None then it's not usable in operations with other quantities. Otherwise it would work as normal, so one just has to do something like:

def __add__(self, other):
    if self.system is None:
        raise ValueError('blah blah blah')
    return super(EMConstant, self).__add__(other)

Thoughts?

@iguananaut - I like your suggestion of checking if system is None and behaving accordingly. The only thing that might be a slight pain is that we make sure we do it for all operators that Quantity overrides.

@iguananaut - do you plan to try implementing this, or should I do it? (I may have a little time this evening)

Yeah, that is a bit of a pain but I don't know if I see a better way :/ I suppose hypothetically I could write a metaclass that does this automatically :)

@iguananaut - one of my ideas was to create a new class called e.g. EMConstantHolder or something like that that doesn't inherit from Quantity, and then has attributes for the different systems, which return an actual EMConstant object. Then EMConstantHolder doesn't have any operations overloaded.

@iguananaut - so just to make sure, are you going to try tackling this, or shall I try later if I get the chance?

I'll go ahead and give it a try and send you a PR.

@iguananaut - thanks!

@iguananaut - just to check, are you still planning on doing this? If we get something that works nicely, I suggest we just go with it and do the release candidate, and there's always time to tweak things slightly. I think there's general agreement that this would be a good way to have the EM constants work.

I'll give it a shot in a few. I'm having second thoughts about calling this a "release candidate". Maybe "beta" would be better since there are still too many loose ends I think. "rc" should be reserved for something closer to done, but where there might still be a few bugs to fix up.

beta sounds good to me!

As I start hacking on this several questions come up regarding the consistency of this implementation. First off, right now e.si returns a Quantity while e.esu and e.gauss return an EMConstant:

>>> e
<Constant: name='Electron charge' value=1.602176565e-19 error=3.5e-27 units='C' reference='CODATA 2010'>
>>> e.si
<Quantity 1.602176565e-19 C>
>>> e.esu
<Constant: name='Electron charge' value=4.80320450571e-10 error=1.049273603e-17 units='Fr' reference='CODATA 2010'>
>>> e.gauss
<Constant: name='Electron charge' value=4.80320450571e-10 error=1.049273603e-17 units='Fr' reference='CODATA 2010'>

So what should it be in these cases? I think maybe e.esu and e.gauss should be pure quantities since their values are being defined in terms of e's SI definition. Or would the idea be to keep those a Constants and require e.esu.cgs and e.esu.si (note: currently the former is broken).

I think e.anything should return a quantity, for consistency. It sounds like that's what you are thinking too?

@astrofrog I think that makes sense -- thanks for working on this @iguananaut!

@iguananaut - just to check, have you had a chance to work on this more? I think we should try and get this settled and back-ported into 2.0.x (and also fix #603 at the same time).

Huh. I wrote a long comment a couple weeks ago when I was looking into this that appears to have disappeared?!

In that comment I wrote that it would make more sense to tackle this after #579 is merged in. I wonder what happened there... maybe I never actually posted it.

@astrofrog I'm seeing now your idea behind an EMConstantHolder class, but ultimately I still think that might make it too much of a special case. For example, something like isinstance(e, Constant) would no longer work. I don't know exactly what the use case is for it, but it still seems odd to suddenly have two separate classes hierarchies that can represent some constant.

Maybe instead we could flip things around a bit--go back to using the Constant class the way it was before we made it a subclass of Quantity. That is, make Constant have the same purpose ConstantDefinition has now. Then define two two types of Constant subclasses (I don't have good names for them as of yet): One for constants that have the same physical units in different systems of units and can be combined with Quantitys--it could have both Quantity and Constant as its subclasses. The other type is for constants whose units depend on the system, such as the EMConstant class we're currently discussing. In fact, for now we could eschew a middle-class and just have EMConstant be a direct subclass of Constant with the appropriate properties defined on it.

Now you FooConstants (still needs a name) that can be freely combined with other quantities, and EMConstants which cannot without first asking for a specific system. Both are instances of Constant. This approach would also be compatible for other systems that are incompatible with SI, such as where G is defined as 1.

Actually, I have yet another idea that's sort of a synthesis of everything I described above. Let's see if it works...

I've almost got this working, but I need to step away for a little bit and focus on getting beta 2 out. We might want to push this to 0.2.1 if my solution ends up being too destructive (in the meantime we can certainly accept what we already have as a temporary solution).

@astrofrog @iguananaut - I think we want this in 0.2 ... it's a pretty major change with a variety of subtle consequences, so it would be better to finalize what we want.

So I guess it depends if we're willing to either do a beta 3 or merge this without a beta.

I've merged #540 to make sure that whatever happens, we at least have a description of the current behavior. I think there's no point in delaying beta 2 for this if you think it'll take a while longer, but if we do include it in 0.2 in the end, we should do a release candidate to make sure it all works fine (and then do the final release very shortly after).

+1 to @astrofrog 's suggestion here of not delaying beta 2 for this, but still getting it in 0.2 . We probably don't generally want major changes between beta 2 and RC, but since we're planning explicitly for this one, its fine.

@iguananaut @astrofrog - now that beta 2 is out, how do we want to proceed? @iguananaut, are you still working on the largish change you mentioned above?

Yeah, sorry, I really went down the rabbit hole with this but I think I've almost got it. Basically where I'm at is that a Constant instance may or may not have a 'system' attached to it. If it has a 'system' attached it has a value and unit in the appropriate system. If it does not have a system attached to it, it will basically use the 'default' system (i.e. the first system it was defined with, usually SI) unless the systems in which the constant is defined use incompatible units. Then trying to combine a constant with other values without specifying a system will result in a TypeError.

There's still also a special EMConstant class whose sole purpose is to raise an exception when converting to .cgs explaining to use .esu or .gauss instead.

How does that all sound?

Seems reasonably sound to me, @iguananaut ... and is that a fresh branch, or does it build off of this.

Specifically, does it preserve the change that you can just do from astropy.constants import G and not have to worry about the system unless it's an EM system where it's not freely convertable? (I think that's in this PR, although I may have gotten it confused with some of the other discussions...)

@iguananaut - that sounds fine to me! Feel free to open a PR against my branch, or open a new PR (either way is fine). I'll +1 @eteq's comment that for most cases, the user should now have to care about these internal aspects.

Right--they should be able to import a constant directly from astropy.constants, but if it happens to have definition in different units that are not convertible then any attempt to use the constant in an expression will result in a TypeError unless a system is specified.

(Note: This means one can manually define a dimensionless version of G, which will cause all instances of G to become dimensionally ambiguous. That may or may not be a problem--I don't know.)

@iguananaut - I'm not sure I understand your last comment, but it will be easier to see once the code is ready (maybe you can provide an example once it is).

@iguananaut - shall we try and finalize this this week? I've rebased this branch in case you want to add your changes to it, but feel free to also open an alternative PR instead. If we think it's going to take a while to get something stable, then I propose we postpone 0.3, since we are so close to releasing 0.2!

I would be strongly in favor of getting this in for 0.2 (i.,e., end of this week). I think it would be very irritating to users to require from astropy.constants import si/from astropy.constants.si import G in 0.2 and then in 0.3 change it to from astropy.constants import G. It'd be fineto delay a complete solution to the EM units case as long as we still get the from astropy.constants import G in now, but I think we don't want to change how people import fundamental constants from 0.2 to 0.3 .

@eteq - note that this PR already implements from astropy.constants import G for fundamental constants, so if @iguananaut thinks it will take a while to implement the improved framework, then we could merge this for 0.2 and make the internal changes for 0.3 which may require a small change in API for the EM constants for users.

There's a failing test and the docs need updating, but in terms of API, I don't think things would change too much with @iguananaut's changes - here's the current behavior:

In [1]: from astropy.constants import G, k_B

In [2]: G
Out[2]: <Constant: name='Gravitational constant' value=6.67384e-11 error=8e-15 units='m3 / (kg s2)' reference='CODATA 2010'>

In [3]: k_B
Out[3]: <Constant: name='Boltzmann constant' value=1.3806488e-23 error=1.3e-29 units='J / (K)' reference='CODATA 2010'>

In [4]: from astropy.constants import e

In [5]: e.si
Out[5]: <Quantity 1.602176565e-19 C>

In [6]: e.gauss
Out[6]: <Constant: name='Electron charge' value=4.803204505713468e-10 error=1.049273603e-17 units='Fr' reference='CODATA 2010'>

In [7]: e.esu
Out[7]: <Constant: name='Electron charge' value=4.803204505713468e-10 error=1.049273603e-17 units='Fr' reference='CODATA 2010'>

In [8]: e.cgs
ERROR: TypeError: Cannot convert EM constants to cgs because there are different systems for E.M constants within the c.g.s system (ESU, Gaussian, etc.). Instead, directly use the constant with the appropriate suffix (e.g. e.esu, e.gauss, etc.). [astropy.constants.constant]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-8-ebc8d90fd7bd> in <module>()
----> 1 e.cgs

/Users/tom/Library/Python/3.2/lib/python/site-packages/astropy-0.3.dev3119-py3.2-macosx-10.6-x86_64.egg/astropy/constants/constant.py in cgs(self)
     53     @property
     54     def cgs(self):
---> 55         raise TypeError("Cannot convert EM constants to cgs because there "
     56                         "are different systems for E.M constants within the "
     57                         "c.g.s system (ESU, Gaussian, etc.). Instead, "

TypeError: Cannot convert EM constants to cgs because there are different systems for E.M constants within the c.g.s system (ESU, Gaussian, etc.). Instead, directly use the constant with the appropriate suffix (e.g. e.esu, e.gauss, etc.).

It's not perfect, in that e.si returns a Quantity and e.esu a Constant, but if it works, maybe it's good enough for 0.2?

My thing is almost done. I'm out sick today, but tomorrow I'll try to push it. It's built on top of this PR and has many similarities from the user's perspective at least.

@iguananaut - ok, thanks! (and get better soon!)

@iguananaut - do you need me to rebase this? (I merged in #740, which I think is causing the conflict)

Uhh...not just yet. I have my own branch of this that I've already rebased several times. I'll probably make a separate PR once it's closer to done (my task for the rest of today I hope).

Ok, great! I'll leave it as is, and a separate PR probably makes sense anyway.

We close this without merging if #745 gets merged, right?

@eteq - yes, that's correct.

And this can be closed! Yay! \o/