Divide by zero when transforming to numerics

[Jason Moore]

In the mechanics package we sometimes end up with symbolic expressions that can give a divide by zero error if we use the expression in a numerical function. Below is a simple example of an expression that could be generated by the mechanics package. It should evaluate to zero (i.e the limit as vx -> zero. evalf says it's zero but if I use lambdify (or probably any code generator) to transform the expression to something that is purely numeric we get divide by zeros. I feel like this must be a common issue when transforming symbolic expressions to numerical functions. Does anyone know how or if this is dealt with? It be nice if we could generate equations of motion in the mechanics package and generate code from them that wouldn't have divide by zero errors, requiring some fix on the numerical side.

In [1]: from sympy import symbols, sqrt, lambdify

In [2]: k, vx, vy, vz = symbols('k vx vy vz')

In [3]: f = k * sqrt(vx ** 2 + vy ** 2 + vz ** 2)

In [4]: dfdvx = f.diff(vx)

In [5]: dfdvx
Out[5]: k*vx/sqrt(vx**2 + vy**2 + vz**2)

In [6]: dfdvx.evalf(subs={k: 1.0, vx: 0.0, vy: 0.0, vz: 0.0})
Out[6]: 0

In [7]: numeric_dfdvx = lambdify((k, vx, vy, vz), dfdvx)

In [8]: numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-8-55e44f8ab1bf> in <module>()
----> 1 numeric_dfdvx(1.0, 0.0, 0.0, 0.0)

/usr/local/lib/python2.7/dist-packages/numpy/__init__.pyc in <lambda>(_Dummy_26, _Dummy_27, _Dummy_28, _Dummy_29)

ZeroDivisionError: float division by zero

In [9]: numeric_dfdvx = lambdify((k, vx, vy, vz), dfdvx, modules='numpy')

In [10]: numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
/usr/local/lib/python2.7/dist-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in double_scalars
  """
Out[10]: nan

In [11]: numeric_dfdvx = lambdify((k, vx, vy, vz), dfdvx, modules='mpmath')

In [12]: numeric_dfdvx(1.0, 0.0, 0.0, 0.0)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-12-55e44f8ab1bf> in <module>()
----> 1 numeric_dfdvx(1.0, 0.0, 0.0, 0.0)

<string> in <lambda>(k, vx, vy, vz)

/usr/local/lib/python2.7/dist-packages/sympy/mpmath/ctx_mp_python.pyc in __rdiv__(s, t)
    206         if t is NotImplemented:
    207             return t
--> 208         return t / s
    209
    210     def __rpow__(s, t):

/usr/local/lib/python2.7/dist-packages/sympy/mpmath/ctx_mp_python.pyc in __div__(self, other)

/usr/local/lib/python2.7/dist-packages/sympy/mpmath/libmp/libmpf.pyc in mpf_div(s, t, prec, rnd)
    928     if not sman or not tman:
    929         if s == fzero:
--> 930             if t == fzero: raise ZeroDivisionError
    931             if t == fnan: return fnan
    932             return fzero

ZeroDivisionError:

In [19]: dfdvx.limit(vx, 0.0)
Out[19]: 0

Jason
moorepants.info
+01 530-601-9791

[Stefan Krastanov]

The poorman's solution that I used when I had this problem was to do
the calculations in a different basis. I do not know how
general/automatic this solution would be.

However, a common reason for bad/no simplification is lack of
appropriate assumptions. Define all you symbols as real/positive and
half of the issues will go away.

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> For more options, visit https://groups.google.com/groups/opt_out.

[Jason Moore]

In this case, real and positive assumptions don't help for the lambdify call. I guess making lambdify assumption aware could be an option.

I think we should define dynamicsymbols in mechanics to have a default real assumption, this could help in simplification.

The different basis is an interesting idea. Do you have an example of that?

Jason
moorepants.info
+01 530-601-9791

[Matthew Rocklin]

I would use a Piecewise to define something like a sync function.  This is an explicit way to do what you want on the SymPy side.

[Jason Moore]

So our module generates tons of .diff() calls on expressions with generally multiplies, divides, powers, adds, subtracts, and the basic three trig functions (cos, sin, tan). These expressions can be really big so we can't identify these things by eye all the time. How would I have a Piecewise map to these subexpressions that would cause a divide by zero?

Jason
moorepants.info
+01 530-601-9791

[Matthew Rocklin]

This is why I limited my statement to "something like a sync function"

[Stefan Krastanov]

> In this case, real and positive assumptions don't help for the lambdify call. I guess making lambdify assumption aware could be an option.

I was actually expecting that one would always use `simplify` before
`lambdify`. This would make a big difference, because with the
appropriate assumptions `simplify` would remove all such removable
singularities[1], which is the technical term for these point.
Unfortunately a quick google search for "numerical removable
singularity" does not show anything promising.

[1]: https://en.wikipedia.org/wiki/Removable_singularity

[Jason Moore]

I'll don't have a concrete example to show at the moment, but I feel like we get these singularity issues when we do use simplification and don't when we don't or when we expand everything out. I'll try to put together an example where this is the case.

Jason
moorepants.info
+01 530-601-9791

[Jason Moore]

Also, why does .evalf give zero in my previous example? Shouldn't that also raise a DivideByZero error?

Jason
moorepants.info
+01 530-601-9791

[Aaron Meurer]

With your example, it matters what order you substitute the variables.
If you substitute vx=0 first, you get 0/sqrt(vy**2 + vz**2) == 0. If
you substitute vy=vz=0, you get vx/sqrt(vx**2), which should simplify
to sqrt(vx) for vx >= 0 (which it is, since vx=0). Can you assume
nonnegative for these variables. Of course, if you substitute all
three at once, you get 0/0.

There's no general solution to this problem, but you aren't working in
a completely general situation. I think the solution you should adopt
will depend on just how much structure you know about your
expressions.

Sometimes you can rewrite expressions so that they don't have 0/0
(like you learn in calculus). I don't see how to do it for this
expression, though.

Aaron Meurer

[Jason Moore]

Thanks Aaron. The main assumption that we can make for our package is that these variables are real numbers, so that will still leave us with divide by zero issues. Maybe once I have a list of expressions that are common to our problems, we can figure out how to deal with that subset.

Jason
moorepants.info
+01 530-601-9791

[Gilbert Gede]

I think a really common example that came up in mechanics was with trig expressions. Here's a simple example (simple relative to actual mechanics expressions):

sin(q1(t))*sin(q2(t))*cos(q1(t))  /  (sin(q1(t))*sin(q2(t)) + sin(q1(t))*cos(q2(t)))

When I substitute in q1(t) = 0, q2(t) = 0, I get nan.

Generally, we don't know enough about the structure to know which order to substitute things in, i.e. of the coordinates (q's), which ones will end up in the denominator or even if you'll have an expression where both variables end up being able to cause a nan:

sin(q1(t))*sin(q2(t))*cos(q1(t))  /  (sin(q1(t))*sin(q2(t)) + sin(q1(t))*cos(q2(t)))  +

sin(q2(t))*sin(q1(t))*cos(q2(t))  /  (sin(q2(t))*sin(q1(t)) + sin(q2(t))*cos(q1(t)))

Some inspection/simplification seems to be required, which isn't practical for the larger expressions we can generate.

Another question we can ask, relevant to this problem in mechanics, is when these nan's that pop up, should the result always be 0? That knowledge might be useful.

-Gilbert

[Jason Moore]

It should only evaluate to zero if the limit of the expression with respect to the evaulated variables goes to zero. 1/x, for example, would not be one that should evaluate to zero. But it would certainly be interesting to know if our operations plus dervitives in mechanics, only ended up with expressions that should evaluate to zero. Maybe we could prove that somehow...

Jason
moorepants.info
+01 530-601-9791

[Aaron Meurer]

In this case, my rewrite idea works. For instance, the result from
trigsimp no longer has the removable singularity

In [10]: trigsimp(a)
Out[10]:
___
╲╱ 2 ⋅sin(q₂(t))⋅cos(q₁(t))
───────────────────────────
⎛ π⎞
2⋅sin⎜q₂(t) + ─⎟
⎝ 4⎠

In [12]: trigsimp(a).subs(S("q1(t)"), 0).subs(S("q2(t)"), 0)
Out[12]: 0

The Fu trigsimp algorithm is pretty customizable if you dig down into
it, so if you can figure out a reliable trigonometric transformation
that fixes this issue, you can apply it without the full power of
trigsimp. Or rewrite your algorithms to avoid them in the first place,
if you can pinpoint the source.

Aaron Meurer