Re: [sympy] Undesirable behavior from cse

[Chris Smith]

On Fri, Sep 7, 2012 at 1:40 AM, Larry Wigton <lbw...@aol.com> wrote:
> Python code using cse from sympy:
>
> from sympy import *
> x=Symbol('x')
> y=Symbol('y')
> eq1 = 5*x**3*y**2 + y**3
> eq2 = 4*x**2*y**3 + y**2
> eq = [eq1,eq2]
> print eq
> (red,rep) = cse(eq)
> print red
> print rep
> eq = [eq2,eq1]
> print eq
> (red,rep) = cse(eq)
> print red
> print rep
>
> ***********************************************
>
> Output from the code:
>
> [5*x**3*y**2 + y**3, 4*x**2*y**3 + y**2]
> [(x0, y**3), (x1, x0**(2/3))]
> [5*x**3*x1 + x0, 4*x**2*x0 + x1]
>
> [4*x**2*y**3 + y**2, 5*x**3*y**2 + y**3]
> [(x0, y**2), (x1, x0**(3/2))]
> [4*x**2*x1 + x0, 5*x**3*x0 + x1]
>
> **********************************************

Current master gives

[5*x**3*y**2 + y**3, 4*x**2*y**3 + y**2]
[(x0, y**2)]
[x0*(5*x**3 + y), x0*(4*x**2*y + 1)]
[4*x**2*y**3 + y**2, 5*x**3*y**2 + y**3]
[(x0, y**2)]
[x0*(4*x**2*y + 1), x0*(5*x**3 + y)]

googling sympy bleeding edge if you would like to get the most
up-to-date (and soon to be released) 0.7.2 version.

[Larry Wigton]

Thank you for telling me about the Bleeding edge version of sympy.

Unfortunately it does not work for case which worked before.

Test case for bleeding edge version of sympy:

**************************************************

pg00 = Symbol('pg00')
pg01 = Symbol('pg01')
eq1 = sympify(pg01**2/pg00**2)
eq2 = sympify(1.0 - 2*pg01/pg00)
eq3 = sympify(pg01*(1.0 - pg01/pg00))
eq = [eq1,eq2,eq3]

print eq
(red,rep) = cse(eq)
print red
print rep

********************************************************

Output from code:

[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(1.0 - pg01/pg00)]
Traceback (most recent call last):
  File "mytest3.py", line 12, in <module>
    (red,rep) = cse(eq)
  File "/acct/lbw9902/gitstuff/sympy/sympy/simplify/cse_main.py", line 313, in cse
    reduced_exprs[j] = update(expr)
  File "/acct/lbw9902/gitstuff/sympy/sympy/simplify/cse_main.py", line 309, in <lambda>
    update = lambda x: x.subs(subtree, sym)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/basic.py", line 852, in subs
    rv = rv._subs(old, new)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/cache.py", line 88, in wrapper
    func_cache_it_cache[k] = r = func(*args, **kw_args)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/basic.py", line 951, in _subs
    rv = fallback(self, old, new)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/basic.py", line 938, in fallback
    arg = arg._subs(old, new, **hints)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/cache.py", line 88, in wrapper
    func_cache_it_cache[k] = r = func(*args, **kw_args)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/basic.py", line 949, in _subs
    rv = self._eval_subs(old, new)
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/mul.py", line 1276, in _eval_subs
    rat.append(ndiv(c_e, old_e))
  File "/acct/lbw9902/gitstuff/sympy/sympy/core/mul.py", line 1186, in ndiv
    if not b.q % a.q or not a.q % b.q:
AttributeError: 'Infinity' object has no attribute 'q'

**************************************************************

Output from sympy 7.1

[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(1.0 - pg01/pg00)]
[(x0, 1/pg00), (x1, pg01*x0)]
[x1**2, -2*x1 + 1.0, pg01*(-x1 + 1.0)]

***********************************************************

So it looks like bleeding edge version of sympy still needs some work.

   Larry Wigton

 

[Ondřej Čertík]

Larry,

This looks like a regression. Thanks for posting it.
If you want to help us debug it, just use "git bisect" to figure out a
patch that broke this.
Then it should be possible to see how to fix this.

Ondrej

[Chris Smith]

>> if not b.q % a.q or not a.q % b.q:
>> AttributeError: 'Infinity' object has no attribute 'q'

Infinity stopped inheriting from Rational a long time ago and this is
one of those lurking areas where it causes problems.

[Ondřej Čertík]

What would be the way to fix this?

Ondrej

[Chris Smith]

I only looked at this very quickly:

I get


>>> from sympy import *

>>> pg00 = Symbol('pg00')
>>> pg01 = Symbol('pg01')
>>> eq1 = sympify(pg01**2/pg00**2)
>>> eq2 = sympify(1.0 - 2*pg01/pg00)
>>> eq3 = sympify(pg01*(1.0 - pg01/pg00))
>>> eq = [eq1,eq2,eq3]
>>> print eq

[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(1.0 - pg01/pg00)]

>>> (red,rep) = cse(eq)
>>> red
[(x0, -pg01/pg00)]
>>> rep
[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(x0 + 1.0)]


after changing ndiv to be


def ndiv(a, b):
"""if b divides a in an extractive way (like 1/4 divides 1/2
but not vice versa, and 2/5 does not divide 1/3) then return
the integer number of times it divides, else return 0.
"""
if a not in (S.Infinity, S.NegativeInfinity) and \
(not b.q % a.q or not a.q % b.q):
return int(a/b)
return 0

But it could be that the better fix is to not have infinity pulled out
-- by breakup?

Hopefully that can get someone started.

/c

[Larry Wigton]

Thanks for the temporary fix to ndiv.

However even with this fix cse in bleeding edge sympy gives:

[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(1.0 - pg01/pg00)]

[(x0, -pg01/pg00)]

[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(x0 + 1.0)]

whereas cse in sympy 7.1 gives:

[pg01**2/pg00**2, 1.0 - 2*pg01/pg00, pg01*(1.0 - pg01/pg00)]

[(x0, 1/pg00), (x1, pg01*x0)]
[x1**2, -2*x1 + 1.0, pg01*(-x1 + 1.0)]

I think sympy 7.1 result is better.

    Larry Wigton

The introduction of the powers 2/3 and 3/2 is not acceptable.

For one thing in Fortran 2/3 will beocme 0 and 3/2 will become 1.

Even if we make 2 and 3 floating point numbers this is not efficient.

Also if we take x0=y**2 and then try to compute y**3 using x0**(1.5)
we get the wrong answer in the case y is negative.

************************************************************

Is there some way to run cse to avoid this undesirable behavior?

Can cse be modified to avoid this undesirable behavior?