# units tests

I just stumbled over a problem again that I reported years ago in sage-devel, but perhaps someone here can help. I am still after the quick and dirty way for using the existing units package in sage to test symbolic equations for consistent units. This is currently not possible because units are treated as conventional variables, meaning that:

units.x.y - units.x.y

0

The only thing that really stops me right now is that I don't know how to get it to return:

units.x.y + units.x.y

units.x.y

units.x.y - units.x.y

units.x.y

Is there an easy way to write a function or create a subclass that would do that while keeping the rest (multiplication, division, exponents etc.) the way it is?

UPDATE: A modification of the first answer by tmonteil does part of the job:

class MYExp(sage.symbolic.expression.Expression):
def __sub__(self,other):
if self.convert() == other.convert():
return self
else:
raise TypeError('Subtraction of non-matching units')

if self.convert() == other.convert():
return self
else:
raise TypeError('Sum of non-matching units')
def __mul__(self,other):
return MYExp(SR, super(MYExp, self).__mul__(other))
def __div__(self,other):
return MYExp(SR, super(MYExp, self).__div__(other))
def __pow__(self,other):
return MYExp(SR, super(MYExp, self).__pow__(other))

U = lambda x: MYExp(SR,x)
for s in units.length.trait_names():
globals()[s] = U(getattr(units.length,s))


Now this works:

type(meter*2)

class '__main__.MYExp'

But this does not:

type(2*meter)

type 'sage.symbolic.expression.Expression'

Similar for division. Is there a way to somehow modify the above to make any expressions involving MYExp return a MYExp, rather than just those expressions that start with MYExp?

edit retag close merge delete

This was explained at the end of my answer : if you want 2*meter lead to a MYExp expression, you have to define a coercion.

( 2015-03-09 16:37:54 +0200 )edit

Adding reverse multiplication and division did the trick (see my proposed solution below). Defining coercions seems outside my area of comfort.

( 2015-03-09 17:43:33 +0200 )edit

Sort by » oldest newest most voted

You can define a class MYExp, that inherits from symbolic expressions (so that you keep all other operations), but owerwrites substraction and addition to avoid cancellations:

class MYExp(sage.symbolic.expression.Expression):
def __sub__(self,other):
if self == other:
return self
else:
def __mul__(self,other):
return MYExp(SR, super(MYExp, self).__mul__(other))
def __div__(self,other):
return MYExp(SR, super(MYExp, self).__div__(other))
U = lambda x: MYExp(SR,x)

watt = U(units.power.watt)
meter = U(units.length.meter)
newton = U(units.force.newton)
second = U(units.time.second)
joule = U(units.energy.joule)


Then, you can do something like:

sage: watt-watt
watt
sage: watt+watt
watt
sage: (watt-watt)*(newton+meter) - (watt-watt)*(newton+meter)
(meter + newton)*watt


But use with caution. You wil have to adapt if you want to use some other operation in your computation:

sage: sqrt(watt) - sqrt(watt)
0
sage: type(sqrt(watt))
<type 'sage.symbolic.expression.Expression'>
sage: watt^2 - watt^2         # __pow__ method is called, not the product by itsef
0
sage: watt*watt - watt*watt
watt^2


Also, some "external" operations use coercion, that should be defined as well:

sage: 2*watt-2*watt
0
sage: watt*2-watt*2
2*watt
sage: type(watt*2)
<class 'sage.all_cmdline.MYExp'>
sage: type(2*watt)
<type 'sage.symbolic.expression.Expression'>


If you want also to deal with such external products, you will have to define a class MYSR at the same level as SR (Symbolic Ring), and define the .parent() method of an element of MYExp to be MYSR, and define the coercion in MYSR. See those pages if you want to deal with such external products :

more

This is awesome, thanks!! Now the way I would like to use it is to save the units of all my symbolic variables in a dictionary (e.g. udict) and substitute them in a given equation. Unfortunately, (var1 - var2).subs(udict) still gives 0. However, if I do udict[var1] - udict[var2], I get the desired units as in your example. Any ideas? UPDATE: Apparently, expression.subs(somethingelse), even if all elements of expression are substituted by a different type, always stays an expression.

( 2013-06-12 04:48:16 +0200 )edit

Now, another possibe point of view (perhaps more pedagogical) could be not to allow adding things with incompatible units (you can use the convert() function to test that: this function transforms your units in the international system standard). For this, you have to modify the .__add__() method in the class i defined above.

def __add__(self,other):
if self.convert() == other.convert():
return self
else:
raise TypeError('Physics do not like such additions')

sage: (joule + second*watt)*newton
joule*newton
sage: joule-joule
joule
sage: joule - second
TypeError: Physics do not like such additions


If you also define .__pow__() method as of above (imitate the product), i guess you have a way to test wether an expression makes sense, and in which units wil be the result.

more

This is exactly what I want to do, and believe me, it is not just pedagogical, but you'd be surprised how many equations in the scientific literature fail this test, only that it is often less obvious than this. Would be good to add something like this to the units package in sage.

( 2013-06-12 05:13:09 +0200 )edit

Works perfectly with the .__pow__() method, but what should I do with sqrt()? Example: joule^(1/2) - joule^(1/2) gives sqrt(joule), but sqrt(joule) - sqrt(joule) gives 0.

( 2013-06-12 11:49:57 +0200 )edit

Concerning your second comment, see the explanation on the previous answer, this is because the .sqrt() method is inherited from Symbolic Expression. Concerning your first comment, if we want to have this included in Sage, this will require much more work (for example, we will have to teach Sage that cos(x) is defined only if x has no physical dimension, define .sqrt() method as above, and so on until the module can eat most of the symbolic expressions, also we will have to write all the related documentation). Do you feel you could work on this ? We can do it together if needed, so we can first work on a common file until we get something good, and then transform it into a patch. Drop me an email if you are interested, and i will open a ticket !

( 2013-06-16 17:37:29 +0200 )edit

An alternative would be to just write a function that dissects an equation into powers, logs, numerators, denominators, terms, checks each term separately and puts it back together. How can I contact you by email?

( 2013-06-19 07:43:16 +0200 )edit

A small edit to tmonteil's first answer allows handling multiplications and divisions in the above framework. Need to define reverse multiplication and division (http://www.rafekettler.com/magicmetho...):

class MYExp(sage.symbolic.expression.Expression):
def __sub__(self,other):
if self.convert() == other.convert():
return self
else:
raise TypeError('Subtraction of non-matching units')

if self.convert() == other.convert():
return self
else:
raise TypeError('Sum of non-matching units')
def __mul__(self,other):
return MYExp(SR, super(MYExp, self).__mul__(other))
def __rmul__(self,other):
return MYExp(SR, super(MYExp, self).__mul__(other))
def __div__(self,other):
return MYExp(SR, super(MYExp, self).__div__(other))
def __rdiv__(self,other):
return MYExp(SR, super(MYExp, self).__div__(other))
def __pow__(self,other):
return MYExp(SR, super(MYExp, self).__pow__(other))

U = lambda x: MYExp(SR,x)
for s in units.length.trait_names():
globals()[s] = U(getattr(units.length,s))


Now we can write:

meter - meter

meter

1*meter - 1*meter

meter

meter^2 - meter^2

meter^2

However, note that external functions can still result in problems, e.g.

sqrt(meter) - sqrt(meter)

0

sin(meter) - sin(meter)

0

Actually, I found a more practical way to automatically check for consistent units. Provided a dictionary containing all variables of eq as keys and each variable multiplied by its units as entries:

def units_check(eq, i0 = 1):
'''
Checks whether all arguments are keys in udict and returns simplified units
If units cancel out on one side while they should not, set i0=2 or something other than 1.
Example:
sage: var('a b c d')
sage: udict = {a: a*units.length.meter^3, b: b*units.length.meter, c: c*units.length.meter, d: d*units.length.meter}
sage: eq = a == b*c*d
sage: units_check(eq)
meter^3 == meter^3
'''
valdict = {}
i = i0
for arg in eq.arguments():
valdict[arg] = i
i = i + 1
udict[arg]
eq_val = eq.subs(valdict)
return (eq.subs(udict).subs(valdict))/eq_val

more

## Stats

Seen: 518 times

Last updated: Mar 09 '15