[SciPy-dev] design of a physical quantities package: seeking comments

[Darren Dale]

On Sunday 03 August 2008 12:54:12 pm Darren Dale wrote:
> On Sunday 03 August 2008 11:35:47 am Jonathan Guyer wrote:
> > On Aug 3, 2008, at 9:22 AM, Darren Dale wrote:
> > >> More, what happens when you have a quantity in Newton-metres (say)
> > >> and
> > >> you ask it to convert to feet? Do you get Newton-feet, or do all the
> > >> occurrences of "feet" in the Newtons get converted?
> > >
> > > They all get converted.
> >
> > I'm with Anne. I'm having a really hard time seeing when I would want
> > this as a default behavior, whereas I frequently want to be able
> > automatically convert BTUs to Joules and atmospheres to Pascals, but I
> > want an exception if somebody tries to give me a pressure when I ask
> > for an energy. Having it possible to do what you suggest seems like it
> > might be useful sometimes, but I don't think it should be the default
> > behavior.
>
> [...]
>
> > Using this scheme:
> >  >>> from fipy import PhysicalField
> >  >>> PhysicalField("12 lyr/cm**3")
> >
> > PhysicalField(12.0,'lyr/cm**3')
> >
> >  >>> PhysicalField("1 lyr") / "7.3 cm**2"
> >
> > PhysicalField(0.13698630136986301,'lyr/cm**2')
> >
> >  >>> PhysicalField("12 lyr/cm**3").inUnitsOf("m**-2")
> >
> > PhysicalField(1.1352876567096959e+23,'1/m**2')
> >
> > so the compound unit that Anne wants is supported and conversion to
> > canonical units happens only when requested.
>
> But in the meantime, you aggregate units like
> "dynes kg ft m^3 s / lbs Km ps^4 watts ohms". Quick, what units should I
> ask for that won't yield an error? It seems like it is more helpful to
> reduce as you go, and later when you want your value expressed with a
> compound unit, you ask for it. Let me try it my way and see what you think,
> and if people are not sold, it should be simple to reimplement so all units
> are aggregated, and reduced only on request.

I think maybe we can have it both ways. I borrowed (and slightly modified) the 
parser enthought put together to interpret unit strings. Perhaps it could be 
modified further such that 'kg*(pc/cm^3)*K/s^4' dictates that pc/cm^3 is a 
compound unit and should not to be automatically decomposed. A method could 
be added to decompose compound units on request, in case that is desired.

I put a mercurial repo up, you can check it out with 
"hg clone \ 
http://dale.chess.cornell.edu/~darren/cgi-bin/hgwebdir.cgi/quantities \ 
quantities"

There is no documentation yet, and I'm a little unhappy with some of the 
organization, but I was really just blitzing through trying to get a demo 
together. It installs in the usual way. There is no documentation, so here is 
a crash course:

In [1]: from quantities import NDQuantity
udunits(3): Already initialized from 
file "/usr/lib/python2.5/site-packages/quan                     
tities/quantities-data/udunits.dat"

In [2]: a=NDQuantity([1,2,3.0],'J')

In [3]: a
Out[3]: NDQuantity([ 1.,  2.,  3.]), kg * m^2 / s^2

In [4]: b=NDQuantity([1,2,3.0],'BTU')

In [5]: b
Out[5]: NDQuantity([ 1055.05585262,  2110.11170524,  3165.16755786]), kg * 
m^2 / s^2

In [6]: c=NDQuantity([1,2,3.0],'m**-2')

In [7]: c
Out[7]: NDQuantity([ 1.,  2.,  3.]), 1 / m^2

In [8]: c.inUnitsOf('parsec/cm^3')
Out[8]: NDQuantity([  3.24077885e-23,   6.48155770e-23,   9.72233655e-23]), 
parsec/cm^3

In [9]: d=c.inUnitsOf('parsec/cm^3')

In [10]: c
Out[10]: NDQuantity([ 1.,  2.,  3.]), 1 / m^2

In [11]: d/a
Out[11]: NDQuantity([  3.24077885e-23,   3.24077885e-23,   3.24077885e-23]), 
s^2 * parsec/cm^3 / kg m^2

In [12]: b/a
Out[12]: NDQuantity([ 1055.05585262,  1055.05585262,  1055.05585262]), 
(dimensionless)


There are some temporary limitations. Enthoughts parser has certain 
expectations, which places some limitations on the units string passed to the 
constructor. udunits, which is used for all unit conversions, has a different 
set of expectations and limitations. This can be improved, I just havent 
gotten around to it yet. The available units are defined in _units.py. 
Quantities cannot yet be created by multiplying a _unit times another number 
or numpy array. It would be better for the unit definitions to be first class 
quantities, then you could do 10*J and get a quantitiy, but there is a 
circular import issue that I need to work through. 

I kept this in for now:

In [3]: a
Out[3]: NDQuantity([ 1.,  2.,  3.]), kg * m^2 / s^2

In [21]: a.units='ft'

In [22]: a
Out[22]: NDQuantity([ 10.76391042,  21.52782083,  32.29173125]), kg * ft^2 / 
s^2

I ran into the concatenate issue that came up a while back, it drops the 
units:

In [26]: numpy.concatenate([a,a])
Out[26]:
array([ 10.76391042,  21.52782083,  32.29173125,  10.76391042,
        21.52782083,  32.29173125])

sqrt doesnt seem to work:

In [28]: numpy.sqrt(a**2+a**2)
Out[28]: NDQuantity([ 15.2224681 ,  30.44493619,  45.66740429]), kg^2 * ft^4 / 
s^4

but this does:

In [30]: (a**2+a**2)**0.5
Out[30]: NDQuantity([ 15.2224681 ,  30.44493619,  45.66740429]), kg * ft^2 / 
s^2

Anyway, feel free to kick the tires and send feedback.


Darren