## Handling units with the quantities module

Posted March 22, 2013 at 10:00 PM | categories: units | tags:

Updated March 23, 2013 at 09:22 AM

The quantities module (https://pypi.python.org/pypi/quantities) is another option for handling units in python. We are going to try the previous example. It does not work, because scipy.optimize.fsolve is not designed to work with units.

import quantities as u import numpy as np from scipy.optimize import fsolve CA0 = 1 * u.mol / u.L CA = 0.01 * u.mol / u.L k = 1.0 / u.s def func(t): return CA - CA0 * np.exp(-k * t) tguess = 4 * u.s print func(tguess) print fsolve(func, tguess)

>>> >>> >>> >>> >>> >>> >>> ... ... >>> >>> >>> -0.00831563888873 mol/L >>> Traceback (most recent call last): File "<stdin>", line 1, in <module> File "c:\Python27\lib\site-packages\scipy\optimize\minpack.py", line 115, in fsolve _check_func('fsolve', 'func', func, x0, args, n, (n,)) File "c:\Python27\lib\site-packages\scipy\optimize\minpack.py", line 13, in _check_func res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) File "<stdin>", line 2, in func File "c:\Python27\lib\site-packages\quantities-0.10.1-py2.7.egg\quantities\quantity.py", line 231, in __array_prepare__ res._dimensionality = p_dict[uf](*objs) File "c:\Python27\lib\site-packages\quantities-0.10.1-py2.7.egg\quantities\dimensionality.py", line 347, in _d_dimensionless raise ValueError("quantity must be dimensionless") ValueError: quantity must be dimensionless

Our function works fine with units, but fsolve does not pass numbers with units back to the function, so this function fails because the exponential function gets an argument with dimensions in it. We can create a new function that solves this problem. We need to “wrap” the function we want to solve to make sure that it uses units, but returns a float number. Then, we put the units back onto the final solved value. Here is how we do that.

import quantities as u import numpy as np from scipy.optimize import fsolve as _fsolve CA0 = 1 * u.mol / u.L CA = 0.01 * u.mol / u.L k = 1.0 / u.s def func(t): return CA - CA0 * np.exp(-k * t) def fsolve(func, t0): 'wrapped fsolve command to work with units' tU = t0 / float(t0) # units on initial guess, normalized def wrapped_func(t): 't will be unitless, so we add unit to it. t * tU has units.' return float(func(t * tU)) sol, = _fsolve(wrapped_func, t0) return sol * tU tguess = 4 * u.s print fsolve(func, tguess)

4.60517018599 s

It is a little tedious to do this, but we might only have to do it once if we store the new fsolve command in a module. You might notice the wrapped function we wrote above only works for one dimensional problems. If there are multiple dimensions, we have to be a little more careful. In the next example, we expand the wrapped function definition to do both one and multidimensional problems. It appears we cannot use numpy.array element-wise multiplication because you cannot mix units in an array. We will use lists instead. When the problem is one-dimensional, the function will take a scalar, but when it is multidimensional it will take a list or array. We will use try/except blocks to handle these two cases. We will assume multidimensional cases, and if that raises an exception because the argument is not a list, we assume it is scalar. Here is the more robust code example.

import quantities as u import numpy as np from scipy.optimize import fsolve as _fsolve def fsolve(func, t0): '''wrapped fsolve command to work with units. We get the units on the function argument, then wrap the function so we can add units to the argument and return floats. Finally we call the original fsolve from scipy. Note: this does not support all of the options to fsolve.''' try: tU = [t / float(t) for t in t0] # units on initial guess, normalized except TypeError: tU = t0 / float(t0) def wrapped_func(t): 't will be unitless, so we add unit to it. t * tU has units.' try: T = [x1 * x2 for x1,x2 in zip(t, tU)] except TypeError: T = t * tU try: return [float(x) for x in func(T)] except TypeError: return float(func(T)) sol = _fsolve(wrapped_func, t0) try: return [x1 * x2 for x1,x2 in zip(sol, tU)] except TypeError: return sol * tU ### Problem 1 CA0 = 1 * u.mol / u.L CA = 0.01 * u.mol / u.L k = 1.0 / u.s def func(t): return CA - CA0 * np.exp(-k * t) tguess = 4 * u.s sol1, = fsolve(func, tguess) print 'sol1 = ',sol1 ### Problem 2 def func2(X): a,b = X return [a**2 - 4*u.kg**2, b**2 - 25*u.J**2] Xguess = [2.2*u.kg, 5.2*u.J] s2a, s2b = fsolve(func2, Xguess) print 's2a = {0}\ns2b = {1}'.format(s2a, s2b)

sol1 = 4.60517018599 s s2a = 2.0 kg s2b = 5.0 J

That is pretty good. There is still room for improvement in the wrapped function, as it does not support all of the options that scipy.optimize.fsolve supports. Here is a draft of a function that does that. We have to return different numbers of arguments depending on the value of full_output. This function works, but I have not fully tested all the options. Here are three examples that work, including one with an argument.

import quantities as u import numpy as np from scipy.optimize import fsolve as _fsolve def fsolve(func, t0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=0.0, factor=100, diag=None): '''wrapped fsolve command to work with units. We get the units on the function argument, then wrap the function so we can add units to the argument and return floats. Finally we call the original fsolve from scipy. ''' try: tU = [t / float(t) for t in t0] # units on initial guess, normalized except TypeError: tU = t0 / float(t0) def wrapped_func(t, *args): 't will be unitless, so we add unit to it. t * tU has units.' try: T = [x1 * x2 for x1,x2 in zip(t, tU)] except TypeError: T = t * tU try: return [float(x) for x in func(T, *args)] except TypeError: return float(func(T)) sol = _fsolve(wrapped_func, t0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag) if full_output: x, infodict, ier, mesg = sol try: x = [x1 * x2 for x1,x2 in zip(x, tU)] except TypeError: x = x * tU return x, infodict, ier, mesg else: try: x = [x1 * x2 for x1,x2 in zip(sol, tU)] except TypeError: x = sol * tU return x ### Problem 1 CA0 = 1 * u.mol / u.L CA = 0.01 * u.mol / u.L k = 1.0 / u.s def func(t): return CA - CA0 * np.exp(-k * t) tguess = 4 * u.s sol1, = fsolve(func, tguess) print 'sol1 = ',sol1 ### Problem 2 def func2(X): a,b = X return [a**2 - 4*u.kg**2, b**2 - 25*u.J**2] Xguess = [2.2*u.kg, 5.2*u.J] sol, infodict, ier, mesg = fsolve(func2, Xguess, full_output=1) s2a, s2b = sol print 's2a = {0}\ns2b = {1}'.format(s2a, s2b) ### Problem 3 - with an arg def func3(a, arg): return a**2 - 4*u.kg**2 + arg**2 Xguess = 1.5 * u.kg arg = 0.0* u.kg sol3, = fsolve(func3, Xguess, args=(arg,)) print'sol3 = ', sol3

sol1 = 4.60517018599 s s2a = 2.0 kg s2b = 5.0 J sol3 = 2.0 kg

The only downside I can see in the quantities module is that it only handle temperature differences, and not absolute temperatures. If you only use absolute temperatures, this would not be a problem I think. But, if you have mixed temperature scales, the quantities module does not convert them on an absolute scale.

import quantities as u T = 20 * u.degC print T.rescale(u.K) print T.rescale(u.degF)

20.0 K 36.0 degF

Nevertheless, this module seems pretty promising, and there are a lot more features than shown here. Some documentation can be found at http://pythonhosted.org/quantities/.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.