moorepants / learn-multibody-dynamics

Interactive computational book on multibody dynamics
https://moorepants.github.io/learn-multibody-dynamics/
Other
121 stars 28 forks source link

Show the value of cse in lambdify now that cse=True is available #121

Closed moorepants closed 1 year ago

moorepants commented 1 year ago

I started drafting this, but hit RecursionErrors in the lambdify compile() call. count_ops() takes a while to compute but only gives 108k operations in expr3. I didn't think 108k operations would be any trouble for lambdify. The operation count has to be reasonably big to see worthwhile timing differences in the resulting functions. I tested Python 3.7 and SymPy 1.9, 1.10 and get the same error. Maybe I'm mistaken about how many operations that lambdify can deal with. I thought it was in the millions.

If we have very large expressions, `lambdify()`` can do use common
subexpressions internally to simplify the numerical code in the generated
function and execute much faster.

.. jupyter-execute::

   a, b, c, x, y, z = sm.symbols('a, b, c, x, y, z')

   expr1 = a*x**2 + a*x**5 + b*sm.sqrt(y) + b*y**4 + c*z + c*(x - y - z)**3
   expr1

.. jupyter-execute::

   expr2 = c*x**2 + a*x**5 + b*sm.sqrt(y) + b*y**4 + c*z + c*(x - y - z)**3
   expr2

.. jupyter-execute::

   expr3 = ((expr1*expr2)**3).expand()
   expr3.count_ops()

.. jupyter-execute::

   eval_expr3_cse = sm.lambdify((a, b, c, x, y, z), expr3, cse=True)

.. jupyter-execute::

   %%timeit
   eval_expr3_cse(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)

.. jupyter-execute::

   eval_expr3 = sm.lambdify((a, b, c, x, y, z), expr3)

.. jupyter-execute::

   %%timeit
   eval_expr3(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
Peter230655 commented 1 year ago

In some of my examples, I have force vectors, where force.count_ops(visual=False) gives 30 mio operations, and no problem with lambdify( [.....], force, cse=True) - except, it may run for 30 min on my machine.

On rare occasions, I did get RecursionsErrors, but it seemed to me I definded some expression stupidly. Frankly I never paid much attention.

Peter230655 commented 1 year ago

I use sympy 1.11.1

moorepants commented 1 year ago

I'm getting:

Python 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:21) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.33.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: run recursion_limit_test.py
1.10.1
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
~/recursion_limit_test.py in <module>
      6 expr2 = c*x**2 + a*x**5 + b*sm.sqrt(y) + b*y**4 + c*z + c*(x - y - z)**3
      7 expr3 = ((expr1*expr2)**3).expand()
----> 8 eval_expr3 = sm.lambdify((a, b, c, x, y, z), expr3)

~/miniconda/envs/sympy1.10/lib/python3.7/site-packages/sympy/utilities/lambdify.py in lambdify(args, expr, modules, printer, use_imps, dummify, cse)
    898     filename = '<lambdifygenerated-%s>' % _lambdify_generated_counter
    899     _lambdify_generated_counter += 1
--> 900     c = compile(funcstr, filename, 'exec')
    901     exec(c, namespace, funclocals)
    902     # mtime has to be None or else linecache.checkcache will remove it

RecursionError: maximum recursion depth exceeded during compilation

In [2]: import numpy as np
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-2-0aa0b027fcb6> in <module>
----> 1 import numpy as np

ModuleNotFoundError: No module named 'numpy'

with this code:

import sympy as sm
print(sm.__version__)

a, b, c, x, y, z = sm.symbols('a, b, c, x, y, z')
expr1 = a*x**2 + a*x**5 + b*sm.sqrt(y) + b*y**4 + c*z + c*(x - y - z)**3
expr2 = c*x**2 + a*x**5 + b*sm.sqrt(y) + b*y**4 + c*z + c*(x - y - z)**3
expr3 = ((expr1*expr2)**3).expand()
eval_expr3 = sm.lambdify((a, b, c, x, y, z), expr3)
moorepants commented 1 year ago

Which was surprising given that expr3.count_ops() returns (slowly) 108k.

Peter230655 commented 1 year ago

I use python version 3.11, no idea whether this matters for lambdify(..)

moorepants commented 1 year ago

We have to use Python 3.8 for this (which was the same as last year) so that it matches our homework cloud system.

Peter230655 commented 1 year ago

In the app I use on my iPad (Carnets Plus), I cannot change the version of python. But even before the latest update, where it was version 3.6 I believe, I never had problems lambdifying expressions with large, 10 mio + operations. I have seen these RecursionErrors, but I always found I had defined some expression 'stupidly'. I changed the definition and it worked. I will pay attention to this in future.

moorepants commented 1 year ago

I wonder if I have something in this particular expression I created that is causing the issue. Maybe I should post to as a SymPy issue.

Peter230655 commented 1 year ago

I tried your example on my 3.11 There is the line: expr3 = ((expr1expr2)*3) .expand() This makes my iPad crash, I do not even get any error message.

If I delete the .expand(), that is if I use expr3 = ((expr1expr2)*3) it works without any problems, it is a harmless term with 41 operations.

As I lambdified expressions with # of operations >> 108 k without any problems, I must assume that this .expand() causes the problem.

Peter230655 commented 1 year ago

My terms are always results of sympy.physics.mechanics, I never tried on 'construed' examples, I would not know how.

Peter230655 commented 1 year ago

If you want, I will send you one of my examples with # of operations >> 100,000 where I have no problem with lambdify(...) >> 100,000

moorepants commented 1 year ago

No need, I have plenty of examples.

For the purposes of this part in the chapter I just wanted a contrived long expression that could be generated in a couple lines of code.

Peter230655 commented 1 year ago

This contrieved code, of around 150,000 operations shows, that lambdifying with cse=True is about 10 times faster than without.

I have seen improvements miuch, much better than this.

import sympy as sm
import time

a, b, c, x, y, z = sm.symbols('a, b, c, x, y, z')
summe = []
for i in range(1000):
    for j in range(10):
        for k in (a, b, c, x, y, z):
            summe.append(k*sm.sin(k)**j)
print('number of operations:', sum([summe[l].count_ops(visual=False) for l in
                                    range(len(summe))]))

summe_lam_cse = sm.lambdify([a, b, c, x, y, z], summe, cse=True)
summe_lam     = sm.lambdify([a, b, c, x, y, z], summe)

start = time.time()
print('cse', sum(summe_lam_cse(1., 2., 3., 4., 5., 6.)))
print('time with cse', time.time() - start)

start = time.time()
print('without cse', sum(summe_lam(1., 2., 3., 4., 5., 6.)))
print('time without cse', time.time() - start)
moorepants commented 1 year ago

Was playing with your example and I hit a MemoryError:

import sympy as sm
a, b, c, x, y, z = sm.symbols('a, b, c, x, y, z')
long_expr = 0

for i in range(500):
    for k in syms:
        expr = (k+i)*sm.sin(k)**(i//10)
        if i % 2 == 0:
            long_expr += expr
        else:
            long_expr *= expr
sm.lambdify(syms, long_expr)
s_push: parser stack overflow
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Cell In[41], line 1
----> 1 sm.lambdify(syms, long_expr)

File ~/miniconda/lib/python3.8/site-packages/sympy/utilities/lambdify.py:888, in lambdify(args, expr, modules, printer, use_imps, dummify, cse)
    886 filename = '<lambdifygenerated-%s>' % _lambdify_generated_counter
    887 _lambdify_generated_counter += 1
--> 888 c = compile(funcstr, filename, 'exec')
    889 exec(c, namespace, funclocals)
    890 # mtime has to be None or else linecache.checkcache will remove it

MemoryError: 

This expression only had 14k operations.

Peter230655 commented 1 year ago

I have no idea! in the many examples I have calculated, at least if they are 3D, the force vector rarely has less than 50,000 operations, often a few hundred thousand, in a few examples up to 60,000,000 operations. I never had any problems with RecursiveError, or stack overflow. ( I did have RecursiveError a handful of times, but it was always related to my defining some term 'stupidly' - I corrected, but paid no attention otherwise)

I wonder what is special with 'contrived' examples.

On Tue 7. Feb 2023 at 21:37 Jason K. Moore @.***> wrote:

Was playing with your example and I hit a MemoryError:

import sympy as sma, b, c, x, y, z = sm.symbols('a, b, c, x, y, z')long_expr = 0 for i in range(500): for k in syms: expr = (k+i)*sm.sin(k)*(i//10) if i % 2 == 0: long_expr += expr else: long_expr = exprsm.lambdify(syms, long_expr)

s_push: parser stack overflow

MemoryError Traceback (most recent call last) Cell In[41], line 1 ----> 1 sm.lambdify(syms, long_expr)

File ~/miniconda/lib/python3.8/site-packages/sympy/utilities/lambdify.py:888, in lambdify(args, expr, modules, printer, use_imps, dummify, cse) 886 filename = '<lambdifygenerated-%s>' % _lambdify_generated_counter 887 _lambdify_generated_counter += 1 --> 888 c = compile(funcstr, filename, 'exec') 889 exec(c, namespace, funclocals) 890 # mtime has to be None or else linecache.checkcache will remove it

MemoryError:

This expression only had 14k operations.

— Reply to this email directly, view it on GitHub https://github.com/moorepants/learn-multibody-dynamics/issues/121#issuecomment-1421411876, or unsubscribe https://github.com/notifications/unsubscribe-auth/AT5MQUSEK34ZXSXYMMC32U3WWKXATANCNFSM6AAAAAAUSYPKU4 . You are receiving this because you commented.Message ID: @.***>

-- Best regards,

Peter Stahlecker