Sympy function support.

[dapperfu]

Trying to use Pint similar to how I used my TI-89 (with units).

Testing the same maths in a Jupyter notebook results in this.

---------------------------------------------------------------------------
SyntaxError                               Traceback (most recent call last)
~/python_gsend/.venv/lib/python3.5/site-packages/sympy/core/sympify.py in sympify(a, locals, convert_xor, strict, rational, evaluate)
    353         a = a.replace('\n', '')
--> 354         expr = parse_expr(a, local_dict=locals, transformations=transformations, evaluate=evaluate)
    355     except (TokenError, SyntaxError) as exc:

~/python_gsend/.venv/lib/python3.5/site-packages/sympy/parsing/sympy_parser.py in parse_expr(s, local_dict, transformations, global_dict, evaluate)
    893 
--> 894     return eval_expr(code, local_dict, global_dict)
    895 

~/python_gsend/.venv/lib/python3.5/site-packages/sympy/parsing/sympy_parser.py in eval_expr(code, local_dict, global_dict)
    806     expr = eval(
--> 807         code, global_dict, local_dict)  # take local objects in preference
    808 

SyntaxError: invalid syntax (<string>, line 1)

During handling of the above exception, another exception occurred:

SympifyError                              Traceback (most recent call last)
<ipython-input-19-83e38e772e91> in <module>()
----> 1 position = integrate(velocity, (t*_s))
      2 position

~/python_gsend/.venv/lib/python3.5/site-packages/sympy/integrals/integrals.py in integrate(*args, **kwargs)
   1289     risch = kwargs.pop('risch', None)
   1290     manual = kwargs.pop('manual', None)
-> 1291     integral = Integral(*args, **kwargs)
   1292 
   1293     if isinstance(integral, Integral):

~/python_gsend/.venv/lib/python3.5/site-packages/sympy/integrals/integrals.py in __new__(cls, function, *symbols, **assumptions)
     73             return function._eval_Integral(*symbols, **assumptions)
     74 
---> 75         obj = AddWithLimits.__new__(cls, function, *symbols, **assumptions)
     76         return obj
     77 

~/python_gsend/.venv/lib/python3.5/site-packages/sympy/concrete/expr_with_limits.py in __new__(cls, function, *symbols, **assumptions)
    354         # This constructor only differs from ExprWithLimits
    355         # in the application of the orientation variable.  Perhaps merge?
--> 356         function = sympify(function)
    357         if hasattr(function, 'func') and function.func is Equality:
    358             lhs = function.lhs

~/python_gsend/.venv/lib/python3.5/site-packages/sympy/core/sympify.py in sympify(a, locals, convert_xor, strict, rational, evaluate)
    354         expr = parse_expr(a, local_dict=locals, transformations=transformations, evaluate=evaluate)
    355     except (TokenError, SyntaxError) as exc:
--> 356         raise SympifyError('could not parse %r' % a, exc)
    357 
    358     return expr

SympifyError: Sympify of expression 'could not parse '1.0 meter / second'' failed, because of exception being raised:
SyntaxError: invalid syntax (<string>, line 1)

[Ricyteach]

You can fix this by extending pint to support the sympify API. See here.

This seems to work ok:

from pint.quantity import _Quantity
from sympy import sympify
_Quantity._sympy_ = lambda s: sympify(f'{s.m}*{s.u:~}')

After which you can do things like this:

>>> import pint
>>> from sympy import *
>>> u = pint.UnitRegistry()
>>> x = 2*u.ft
>>> N(x)
2.0*ft

[hgrecco]

Shall we implement it by default?

[Ricyteach]

I would personally love that. One problem might be that the implementation has to make a choice regarding whether the pint units will be represented in sympy using long from, or using abbreviated units.

Details:

In my example code I used the ~P format specification. The P spec causes the result to be displayed using the pretty print formatting, which is needed to get them to work in sympify, so there is no decision being made there. The ~ is what causes the abbreviated units, so that's optional.

I always use abbreviated units, and it would be annoying to me to have this implemented using the long form version. But others might have a different opinion. Maybe there can be some kind of setting somewhere to tell a unit system to default to abbreviated units always...? That's probably another can of worms.

[hgrecco]

Regarding implementing the simpy API, I think is the way to go because it does no collide with anything else and if simpy is not installed nothing will be broken. (A simpify method will need to be added to compat)

Regarding the format, I think we can make use of the default_format attribute of the UnitRegistry/Quantity/Unit (It might need some improvements)

[Ricyteach]

Aha, I'm not familiar with default_format- must have missed that part of the tutorial. Looks like it would definitely come close to solving the problem.

[Ricyteach]

Using the default_format setting seems to work pretty flawlessly, whether using pretty print or not:

>>> from pint.quantity import _Quantity
>>> from sympy import sympify
>>> _Quantity._sympy_ = lambda s: sympify(f'{s.m}*{s.u}')
>>> import pint
>>> from sympy import *
>>> u = pint.UnitRegistry(system='US')
>>> u.default_format = '~P'
>>> N(2*u.ft**2)
2.0*ft²
>>> u = pint.UnitRegistry(system='US')
>>> N(2*u.ft**2)
2.0*foot**2

[ricopicone]

Was this ever implemented?

[kompre]

Not that I know about.

I delved a bit into the problem, and it is more nuanced than the above method would suggest.

In my experience the method described above will transform the pint unit into a sympy symbol: while this kinda looks ok and works, sympy expression will be not able to understand that some quantity can be simplified or are compatible with each others.

For example imagine you want to add two compatible quantities, but expressed in different units:

1*u.m + 1*u.cm 

sympy will not be able to know the symbol m is worth 100 times the symbol cm, thus you'll end up with a silly looking result. The more complex the expression becomes, the more broken it will looks.

If you want to let sympy to properly recognize a quantity, you'll need to convert pint units into sympy units (sympy.core.physics). Personally I think there is some silliness in that module, because I find particular difficult to define a derived unit, but I'm not a developer so I bet they know better.

I'll paste a piece of code that I use to covert pint units to sympy units. Be aware that this works quite well with common units that are also defined in sympy.core.physics, but doesn't handle well uncommon units (e.g. kgf) or units that may have a different name in sympy than the one used by pint (e.g. gravity).

You'll also need to understand that sympy will not automatically simplify an expression with different units, but you'll need to use the converto_to() function from sympy.core.physics to get the output in the desired units.

Lastly I've tried the code mostly with force (kN), length (m), and pressure (MPa) dimensions since these are in my field of work.

import pint
import sympy.physics.units as sympy_units
from sympy.physics.units.util import convert_to

from sympy import nsimplify, sympify

unitregistry = pint.UnitRegistry()

def pint_to_sympy(quantity: unitregistry.Quantity):
    """convert pint quantity to sympy quantity

    Args:
        quantity (UnitRegistry.Quantity): a quantity defined with the pint module

    """
    # divide and extract the magnitude from the units: it will generate a two elements tuple, where the first item will be the magnitude and the second ona a tuple of tuples; each nested tuple is composed by two elements, the unit proper and the exponent to which is elevated; the tuples are supposed to be multiplied together.

    magnitude, units = (1*quantity).to_tuple() # quantity is multiplied by 1 so that it is converted to pint.Quantity if pint.Unit is passed instead

    # for each unit (i.e. tuple), check if it exist in the sympy.physics.units module
    for u in units:
        fullname = u[0]
        shortname = f'{pint.Unit(fullname):~}'
        exponent = sympify(u[1])

        # add a new unit if it doesn't exist
        if not hasattr(sympy_units, fullname):
            if [True for x in unitregistry.parse_unit_name(fullname) if not x[0] == '']:
                is_prefixed = True
            else:
                is_prefixed = False

            setattr(sympy_units, fullname, sympy_units.Quantity(fullname, abbrev = shortname, is_prefixed=is_prefixed))  # create thwith full namee new sympy unit 
            setattr(sympy_units, shortname, getattr(sympy_units, fullname)) # create the alias with the shortname

            # set the global scale factor relative to base units (base units are assumed to be in sympy)
            # _magnitude, _units = (1 * pint.Unit(fullname)).to_base_units().to_tuple()
            # _reference = sympify(1)
            # for _u in _units:
            #     _reference *= getattr(sympy_units, _u[0])**nsimplify(_u[1])

            # getattr(sympy_units, fullname).set_global_relative_scale_factor(_magnitude, _reference)

        # multiply magnitude for the sympy units (create a sympy.core.Mul object)
        magnitude *= getattr(sympy_units, fullname)**(exponent)

    return magnitude

pint.Quantity._sympy_ = lambda x: pint_to_sympy(x)
pint.Unit._sympy_ = lambda x: pint_to_sympy(1*x)

if __name__ == "__main__":
    u = unitregistry

    F = 5000 * u.daN  # this unit is not present in sympy.core.physics
    A = 2 * u.m
    B = 300 * u.cm

    # pint unit get converted to sympy units
    print(sympify(F))
    print(sympify(A))
    print(sympify(B))

    # sympy will not automatically simplify different units
    print(sympify(F / (A * B)))  # this is a pressure

    # you need to use convert to

    print(convert_to(sympify(F), u.kN))

    print(convert_to(sympify(F / (A * B)), u.MPa))

    print(convert_to(sympify(F / (A * B)), u.kPa))