Wrong calculation of negative powers of variables (Sorry for not describing the error in detail).
itumekanik opened this issue · 10 comments
from numexpr import evaluate
x = 2
print(evaluate("x**(2.0)"))
print(evaluate("x**(-2.0)"))
print(x**(-2))
OUTPUT:
4
0 # WRONG!!!!!
0.25
Is integer truncation playing a role here?
Note that this is only an issue when x is an integer:
>>> from numexpr import evaluate
>>> x = 2
>>> print(evaluate("x**(2.0)"))
4
>>> print(evaluate("x**(-2.0)"))
0
>>> x = 2.0
>>> print(evaluate("x**(2.0)"))
4.0
>>> print(evaluate("x**(-2.0)"))
0.25Hi,
Then why not make "x" float automatically?
Then why not make "x" float automatically?
My guess would be that nobody noticed it (AFAIK numexpr is rarely used to do scalar computations). But I just noticed the problem with negative powers is larger than I thought:
>>> import numexpr as ne
>>> a = np.array([1, 2, 3])
>>> ne.evaluate("a ** -2") # <-- this is correct (mimics numpy)
ValueError: Integers to negative integer powers are not allowed.
>>> ne.evaluate("a ** -2.0") # <-- this is not correct
array([1, 0, 0], dtype=int32)'While with numpy, this gives:
>>> a ** -2
ValueError: Integers to negative integer powers are not allowed
>>> a ** -2.0
array([ 1. , 0.25 , 0.11111111])If that is your case, you may be able to work around the issue by changing your expression (I guess this is what numexpr should do behind the scene -- but I am unsure it would preserve numeric accuracy for high values):
>>> ne.evaluate("1 / (a ** 2.0)")
array([ 1. , 0.25 , 0.11111111])
>>> ne.evaluate("1 / (a ** 2)")
array([ 1. , 0.25 , 0.11111111])Or by casting to float before calling numexpr:
>>> a = a.astype(float)
>>> ne.evaluate("a ** -2.0")
array([ 1. , 0.25 , 0.11111111])
>>> ne.evaluate("a ** -2")
array([ 1. , 0.25 , 0.11111111])I'll be back from vacation on the 22nd, if someone wants to remind me of this issue then I can take a deeper look at it to see if it's resolvable or not.
I think it shouldn't be too hard to solve.
I don't have time to test that, but it seems to me that it would be enough to add a kind argument to the line at:
numexpr/numexpr/expressions.py
Line 315 in 62cc051
Or maybe replace that line with:
r = truediv_op(ConstantNode(1), r)While looking at the code, I think line 319 is also problematic for integer input (i.e. a ** -1.0 is wrong when a is an integer array) and should either use truediv_op or specify kind.
Here is an additional test case for that:
>>> import numexpr as ne
>>> a = np.array([1, 2, 3])
>>> a ** -1.0
array([ 1. , 0.5 , 0.33333333])
>>> ne.evaluate("a ** -1.0")
array([1, 0, 0], dtype=int32)Hi @gdementen ,
I tested your change proposal and replaced line 315 with the following statement
r = truediv_op(ConstantNode(1), r)
This seems to work for my case.
Still I guess it would require some more testing.
Following is a test script of my case for the calculation of two different matrices, one with numexpr the other with numpy. It seems to work fine after the change you mentioned in line 315 of the expressions.py.
import numpy as np
def _lambdifygenerated0(sym_en_3, sym_en_2, sym_en_0, sym_en_1):
from numexpr import evaluate
x0 = sym_en_0 - sym_en_1
x1 = evaluate("x0**2", truediv=True)
x2 = sym_en_2 - sym_en_3
x3 = evaluate("x2**2", truediv=True)
x4 = evaluate("x1 + x3", truediv=True)
x5 = evaluate("x4**1.5", truediv=True)
# x6 = evaluate("5040000.0/x5", truediv=True)
x6 = 5040000.0/x5
x7 = evaluate("x1*x6", truediv=True)
# x8 = evaluate("sqrt(x4)", truediv=True)
x8 = np.sqrt(x4)
x9 = evaluate("sym_en_0**2 - 2*sym_en_0*sym_en_1 + sym_en_1**2 + sym_en_2**2 - 2*sym_en_2*sym_en_3 + sym_en_3**2", truediv=True)
x10 = evaluate("x9**(-2.0)", truediv=True)
x11 = evaluate("x4**1.0", truediv=True)
x12 = evaluate("x9**(-2.5)", truediv=True)
x13 = evaluate("5443200.0*x10*x8 - 10886400.0*x11*x12 + 7257600.0*x5*x9**(-3.0)", truediv=True)
# x13 = 5443200.0*x10*x8 - 10886400.0*x11*x12 + 7257600.0*x5*x9**(-3.0)
x14 = evaluate("x11**(-1.0)", truediv=True)
x15 = evaluate("x13*x14*x3 + x7", truediv=True)
x16 = evaluate("x0*x2", truediv=True)
x17 = evaluate("x16*x6", truediv=True)
x18 = evaluate("x14*x16", truediv=True)
x19 = evaluate("-x13*x18 + x17", truediv=True)
x20 = evaluate("x9**(-1.5)", truediv=True)
x21 = evaluate("x20*x8", truediv=True)
x22 = evaluate("x10*x11", truediv=True)
x23 = evaluate("3628800.0*x12*x5", truediv=True)
x24 = evaluate("3628800.0*x21 - 6350400.0*x22 + x23", truediv=True)
x25 = evaluate("x8**(-1.0)", truediv=True)
x26 = evaluate("x2*x25", truediv=True)
x27 = evaluate("-x24*x26", truediv=True)
x28 = evaluate("-x13", truediv=True)
x29 = evaluate("x14*x28*x3 - x7", truediv=True)
x30 = evaluate("-x17 - x18*x28", truediv=True)
x31 = evaluate("1814400.0*x21 - 4536000.0*x22 + x23", truediv=True)
x32 = evaluate("-x26*x31", truediv=True)
x33 = evaluate("x3*x6", truediv=True)
x34 = evaluate("x1*x14", truediv=True)
x35 = evaluate("x13*x34 + x33", truediv=True)
x36 = evaluate("x0*x25", truediv=True)
x37 = evaluate("x24*x36", truediv=True)
x38 = evaluate("x28*x34 - x33", truediv=True)
x39 = evaluate("x31*x36", truediv=True)
x40 = evaluate("x8*x9**(-1.0)", truediv=True)
x41 = evaluate("x11*x20", truediv=True)
x42 = evaluate("1814400.0*x10*x5", truediv=True)
x43 = evaluate("-x24", truediv=True)
x44 = evaluate("-x26*x43", truediv=True)
x45 = evaluate("x36*x43", truediv=True)
x46 = evaluate("1209600.0*x40 - 2721600.0*x41 + x42", truediv=True)
x47 = evaluate("-x31", truediv=True)
x48 = evaluate("-x26*x47", truediv=True)
x49 = evaluate("x36*x47", truediv=True)
return np.array([[x15, x19, x27, x29, x30, x32],
[x19, x35, x37, x30, x38, x39],
[x27, x37, 2419200.0*x40 - 3628800.0*x41 + x42, x44, x45, x46],
[x29, x30, x44, x15, x19, x48],
[x30, x38, x45, x19, x35, x49],
[x32, x39, x46, x48, x49, 604800.0*x40 - 1814400.0*x41 + x42]])
def _lambdifygenerated1(sym_en_0, sym_en_1, sym_en_2, sym_en_3):
x0 = sym_en_0 - sym_en_1
x1 = x0**2
x2 = sym_en_2 - sym_en_3
x3 = x2**2
x4 = x1 + x3
x5 = x4**1.5
x6 = 5040000.0/x5
x7 = x1*x6
x8 = np.sqrt(x4)
x9 = sym_en_0**2 - 2*sym_en_0*sym_en_1 + sym_en_1**2 + sym_en_2**2 - 2*sym_en_2*sym_en_3 + sym_en_3**2
x10 = x9**(-2.0)
x11 = x4**1.0
x12 = x9**(-2.5)
x13 = 5443200.0*x10*x8 - 10886400.0*x11*x12 + 7257600.0*x5*x9**(-3.0)
x14 = x11**(-1.0)
x15 = x13*x14*x3 + x7
x16 = x0*x2
x17 = x16*x6
x18 = x14*x16
x19 = -x13*x18 + x17
x20 = x9**(-1.5)
x21 = x20*x8
x22 = x10*x11
x23 = 3628800.0*x12*x5
x24 = 3628800.0*x21 - 6350400.0*x22 + x23
x25 = x8**(-1.0)
x26 = x2*x25
x27 = -x24*x26
x28 = -x13
x29 = x14*x28*x3 - x7
x30 = -x17 - x18*x28
x31 = 1814400.0*x21 - 4536000.0*x22 + x23
x32 = -x26*x31
x33 = x3*x6
x34 = x1*x14
x35 = x13*x34 + x33
x36 = x0*x25
x37 = x24*x36
x38 = x28*x34 - x33
x39 = x31*x36
x40 = x8*x9**(-1.0)
x41 = x11*x20
x42 = 1814400.0*x10*x5
x43 = -x24
x44 = -x26*x43
x45 = x36*x43
x46 = 1209600.0*x40 - 2721600.0*x41 + x42
x47 = -x31
x48 = -x26*x47
x49 = x36*x47
return np.array([[x15, x19, x27, x29, x30, x32], [x19, x35, x37, x30, x38, x39], [x27, x37, 2419200.0*x40 - 3628800.0*x41 + x42, x44, x45, x46], [x29, x30, x44, x15, x19, x48], [x30, x38, x45, x19, x35, x49], [x32, x39, x46, x48, x49, 604800.0*x40 - 1814400.0*x41 + x42]])
s0, s1, s2, s3 = 1, 2, 3, 4
a = _lambdifygenerated0(sym_en_0=s0, sym_en_1=s1, sym_en_2=s2, sym_en_3=s3)
b = _lambdifygenerated1(sym_en_0=s0, sym_en_1=s1, sym_en_2=s2, sym_en_3=s3)
print(a-b)
from numexpr import evaluate
x = 2
print(evaluate("x**(2.0)"))
print(evaluate("x**(-2.0)"))
print(x**(-2))Should be fixed, thanks to the contributors.