First time here? Check out the FAQ!

Ask Your Question
1

substitute() raises exception where python does not

asked 8 years ago

stan gravatar image

updated 8 years ago

In my usual workflow, I derive symbolic equations and then substitute their arguments by values stored in dictionaries to obtain numerical solutions. Now I found out that the substitute() command raises division by zero exception where standard python does not. This is a bit of a problem, as these exceptions prevent my code from evaluating all data sets once such an exception occurs and I don't want to wrap each line of code in a try-except construct. Does anyone know how the exception could be avoided? Here is an example:

var('y a x')
eq_y = y == a/x
vdict = {}
vdict[a] = 1.
vdict[x] = 0.
vdict[a]/vdict[x]

+infinity

However, if I substitute vdict into eq_y, I get an exception:

eq_y.subs(vdict)

Traceback (click to the left of this block for traceback)
...
ValueError: power::eval(): division by zero
Preview: (hide)

Comments

You forgot to declare 'y' as a variable...

vdelecroix gravatar imagevdelecroix ( 8 years ago )

3 Answers

Sort by » oldest newest most voted
0

answered 8 years ago

kcrisman gravatar image

updated 8 years ago

I see. You are using a symbolic equation in the latter case, which will definitely barf at dividing by zero, but the things you are putting in are "real literals".

sage: type(a)
<type 'sage.symbolic.expression.Expression'>
sage: type(vdict[a])
<type 'sage.rings.real_mpfr.RealLiteral'>

As such, they obey 'usual' division rules for floating-point things like

sage: vdict[x]/vdict[x]
NaN

So if you don't want mathematically correct behavior, I suggest not using a symbolic variable!

sage: eq_y = y == z
sage: def mysubs(a=.1,x=0.,eq=eq_y):
....:     t = a/x
....:     return eq.subs(z=t)
....: 
sage: mysubs()
y == +infinity
sage: mysubs(.0,.0)
y == NaN

Indeed, we get

sage: vdict = {}
sage: vdict['a'] = 0.
sage: vdict['x'] = 0.
sage: mysubs(**vdict)
y == NaN
Preview: (hide)
link

Comments

Thanks for this, but this involves re-writing the original symbolic equation inside the function, so it defeats my purpose of using a symbolic equation in the first place. Anyway, the explanations given above will help me formulate the real problem, which is how to catch such exceptions in a code while continuing with the evaluation. I'll post it shortly.

stan gravatar imagestan ( 8 years ago )
0

answered 8 years ago

nbruin gravatar image

You can use fast_callable to compile a symbolic expression into something that evaluates more in the way of a python expression:

sage: fc=fast_callable(eq_y.rhs(),vars=[x,a])
sage: fc(0.0,1.0)
+infinity
Preview: (hide)
link
0

answered 8 years ago

vdelecroix gravatar image

Hello,

The problem comes from the fact that Sage floating point number behave differently than other numbers with respect to division by zero. Namely

sage: 1 / 0
Traceback (most recent call last):
...
ZeroDivisionError: rational division by zero
sage: 1. / 0.
+infinity

I opened a discussion. But I think that the second example should raise the same error as the first one. See this sage-devel thread.

Preview: (hide)
link

Comments

Incorrect, at least not if we want to stay sort of close to IEEE specs for floating point stuff. See e.g. http://stackoverflow.com/questions/14...

kcrisman gravatar imagekcrisman ( 8 years ago )

Of course, http://www.walkingrandomly.com/?p=5152 points this out: sage: float(1.0)/float(0.0) so I don't know what we should do... I'd go with IEEE on this, since I assume the MPFR people know what they are doing - ?

kcrisman gravatar imagekcrisman ( 8 years ago )

It is also Python choice

$ python -c "1. / 0."
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ZeroDivisionError: float division by zero
vdelecroix gravatar imagevdelecroix ( 8 years ago )

OK, thanks, I didn't realise that this is an inconsistency in sage (not throwing an exception for 1./0.). I agree that it should throw an exception, as it is not clear whether it should be +Infinity or -Infinity, since 0 is neither positive nor negative.

stan gravatar imagestan ( 8 years ago )

@vdelecroix - your example is what I was trying to show, float(1.0)gives native Python floats as in your example. Not sure why Python doesn't conform to IEEE, despite reading several sources. I believe Numpy behaves "correctly" in this case, see http://www.walkingrandomly.com/?p=515...

kcrisman gravatar imagekcrisman ( 8 years ago )

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 8 years ago

Seen: 626 times

Last updated: Mar 09 '16