Symbolic simplification of large trig expressions
Vinzenz Bargsten
I'd like to symbolically simplify trigonometric expressions, which may
become quite large.
In another thread Aaron Meurer pointed out that
expand(expand(a.rewrite(exp)).rewrite(cos)) may be a solution.
The result is good for the small example expression (takes reasonable
time). But for larger expressions the expand after rewrite(exp) won't
finish or take too much memory. I tried to disable caching, to split the
expression with e.g. Add.make_args and to disable some of expand's hints
(also force=True or using expand_mul). That was rather not successful yet.
Here is example code:
http://www.ovgu.de/bargsten/download/simptest.tar.bz2 (8KB)
What would you recommend?
Note that I have problems to run this code with my python3 sympy
installation (immediately recursion limit exceeded). This is not a
problem with python 2.7 .
Regards,
Vinzenz
Aaron Meurer
Unfortunately, as you've discovered, expand() is not exactly the
fastest part of SymPy. The real solution here would be to implement a
better trigsimp algorithm. The way you're doing it will result in
smaller expressions, but it's inefficient, as each replacement will
replace single terms with the sum of terms, requiring expansion to
simplify.
Disabling caching might help with the memory, but it could make things
go so slow that it's not worth it (especially expand, which uses it
heavily).
Looking at your expressions, I think the main simplifications that are
happening are these:
a*sin(x)**2 + a*cos(x)**2 => a
a*sin(x)*cos(x) => a/2*sin(2*x)
a*cos(x)**2 - a*sin(x)**2 => a*cos(2*x)
Would you agree?
Note that I've written these in a form that makes them useful as
transformations.
I think you may want to try writing a custom simplification routine
that uses these. It would go something like this:
First, expand the original expression and get a list of the monomails.
You can use poly() to do this. For example:
In [100]: a = expressions.smallExpr
If you don't know before hand which trig functions are in your
expression, you can use something like this:
In [101]: from sympy.functions.elementary.trigonometric import
TrigonometricFunction
In [102]: a.atoms(TrigonometricFunction)
Out[102]: set([cos(q4), sin(q1), sin(q4), cos(q1), cos(q2), sin(q2),
sin(q3), cos(q3)])
In [116]: p = Poly(a, sin(q1), sin(q2), sin(q3), sin(q4), cos(q1),
cos(q2), cos(q3), cos(q4), sin(2*q1), sin(2*q2), sin(2*q3), sin(2*q4),
cos(2*q1), cos(2*q2), cos(2*q3), cos(2*q4))
In [117]: p.monoms()
Out[117]:
[(4, 4, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 3, 2, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 3, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 2, 3, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 2, 2, 1, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 2, 2, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 2, 1, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 2, 1, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0),
(4, 2, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),...
Note that I only send the generators to poly() that are needed. dq1,
etc. are treated as constants. I include the double argument terms
because these will be needed after the transformations.
This list means that you have a terms with the given monomials. For
example, (4, 4, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) means that
there is a term sin(q1)**4*sin(q2)**4*sin(q3)*cos(q4). You can obtain
the coefficient of this term by using p.coeffs() (they are in the same
order as the monomials, you can put them together with zip).
The idea would be to work on the monomials, which is much more space
efficient than working on the expressions directly. The above rules
would then be something like (sorry, I've switched from numbering
starting at 1 to numbering starting at 0):
== a*sin(x)**2 + a*cos(x)**2 => a ==
- If (2 + s0, s1, s2, s3, c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0,
dc1, dc2, dc3) and (s0, s1, s2, s3, 2 + c0, c1, c2, c3, ds0, ds1, ds2,
ds3, dc0, dc1, dc2, dc3) are both monomails with the same coefficient
a, then they can be replaced with a single monomial (s0, s1, s2, s3,
c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2, dc3) with
coefficient a.
- If (s0, 2 + s1, s2, s3, c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0,
dc1, dc2, dc3) and (s0, s1, s2, s3, c0, 2 + c1, c2, c3, ds0, ds1, ds2,
ds3, dc0, dc1, dc2, dc3) are both monomails with the same coefficient
a, then they can be replaced with a single monomial (s0, s1, s2, s3,
c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2, dc3) with
coefficient a.
...
These rules can probably be generalized in to a single rule using some
more advanced logic. Note how I allow for things like
sin(q1)**4*cos(q2)**2 + cos(q1)**2*sin(q1)**2*cos(q2)**2 =>
sin(q1)**2*cos(q2)**2. 2 + s0 would be implemented as "if s0 >= 2:".
== a*sin(x)*cos(x) => a/2*sin(2*x) ==
- If (n, s1, s2, s3, n, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2,
dc3) is a monomial with coefficient a, it can be replaced with (0, s1,
s2, s3, 0, c1, c2, c3, ds0 + n, ds1, ds2, ds3, dc0, dc1, dc2, dc3)
with coefficient a/2**n.
- If (s0, n, s2, s3, c0, n, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2,
dc3) is a monomial with coefficient a, it can be replaced with (0, s1,
s2, s3, 0, c1, c2, c3, ds0, ds1 + n, ds2, ds3, dc0, dc1, dc2, dc3)
with coefficient a/2**n.
...
Note that you have to use a/2**n (not a/2!).
== a*cos(x)**2 - sin(x)**2 => a*cos(2*x) ==
- If (s0, s1, s2, s3, 2 + c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0,
dc1, dc2, dc3) is a monomial with coefficient a (2 + s0, s1, s2, s3,
c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0, dc1, dc2, dc3) is a monomial
with coefficient -a, then they can be replaced with a single monomial
(s0, s1, s2, s3, c0, c1, c2, c3, ds0, ds1, ds2, ds3, dc0 + 1, dc1,
dc2, dc3) with coefficient a.
...
The algorithm would then be to iterate these rules until no changes
can be made. The resulting expression is the fully simplified one.
You can probably be smart about the order that the rules are applied
in order to minimize the number of passes.
Be aware that I may have made typos or small logical errors above, but
you get the idea. This algorithm could probably be generalized and
put into trigsimp(), but even a slightly more custom version that fits
your needs should not be too hard to write, and would be much faster
(the only thing that will bog it down memory wise will be the
expansion of the original expression).
It also could be generalized to handle 2*n angle arguments by
repeating the double angle rules over the functions that already have
double angles, if you find that is necessary. And as you can imagine,
it can also be generalized in a bunch of other directions as well.
This logic could probably also be applied by using some fancy replace
calls, but the monomial version will be more precise, less bug prone
(because replace logic doesn't always work right for things like x**6
=> x**2*x**4), will be faster, and will use much less memory.
Sorry if this is more of a description of a solution rather than a
solution itself. I don't really feel like implementing this right
now, though perhaps someone else will (potential GSoC students: this
would be a great way to fulfill your patch requirement).
Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
>
Simon
Vinzenz Bargsten
> Wow, that's a big expression!
>
> Unfortunately, as you've discovered, expand() is not exactly the
> fastest part of SymPy. The real solution here would be to implement a
> better trigsimp algorithm. The way you're doing it will result in
> smaller expressions, but it's inefficient, as each replacement will
> replace single terms with the sum of terms, requiring expansion to
> simplify.
>
> Disabling caching might help with the memory, but it could make things
> go so slow that it's not worth it (especially expand, which uses it
> heavily).
> Looking at your expressions, I think the main simplifications that are
> happening are these:
>
> a*sin(x)**2 + a*cos(x)**2 => a
> a*sin(x)*cos(x) => a/2*sin(2*x)
> a*cos(x)**2 - a*sin(x)**2 => a*cos(2*x)
>
> Would you agree?
to implement it.
If this is working, I can finally get rid of Mathematica and the complete
modeling procedure as well as the required model transformations
for a 7-axes robot (Kuka LWR) are carried out by sympy.
Regards,
Vinzenz
Ondřej Čertík
On Mon, Feb 27, 2012 at 4:27 PM, Vinzenz Bargsten <vbar...@freenet.de> wrote:
[...]
> Thank you very much. I think I get the point, but I will need some time to
> implement it.
> If this is working, I can finally get rid of Mathematica and the complete
> modeling procedure as well as the required model transformations
> for a 7-axes robot (Kuka LWR) are carried out by sympy.
That'd be really cool. Yes, the trigonometric simplification is pretty hard.
Besides that, how is SymPy doing otherwise in terms of speed compared
to Mathematica for your problem?
Ondrej
Vinzenz Bargsten
Thanks for the hint. It seems that these substitutions have different
names or no name is given at all in robotics related literature.....
Regards,
Vinzenz
Vinzenz Bargsten
your question is not easy to answer. It depends on the actual task and
how you implement it / the algorithm. The simplification problem is an
example: Using the expand()-replace()-solution is probably slower
compared to Mathematica's Simplify[], but using the monomial-based
solution may be faster.
On average, I would say sympy is at least as fast as Mathematica for my
application.
The advantages of sympy in my eyes compared to Mathematica are that
sympy is deterministic and you get error messages and the messages point
to the problem. Besides the math functions (in many cases superior or
unique) it has many convenient functions (such as the printing functions
to generate c-code of your expressions) and one can use all other python
functions. -- and sympy is free ;)
Regards,
Vinzenz
Aaron Meurer
Well, in general, you should be able to just call simplify() (or
trigsimp()), and it should just work. So simplification (especially
trig simplification) is one area where we can definitely use
improvement.
>
> The advantages of sympy in my eyes compared to Mathematica are that sympy is
> deterministic and you get error messages and the messages point to the
> problem.
This is the first time I've heard this particular reason for SymPy
being better. How is Mathematica non-deterministic? I've never used
it, except for WolframAlpha, because, as you noted, it's not free, and
I don't own a copy.
I suppose the better error messages are a result of Python including
the traceback?
> Besides the math functions (in many cases superior or unique)
Cool. What functions are unique? I only know of one that I'm pretty
sure is unique (or at least I couldn't find it in their docs), but
it's part of the polys module. I guess maybe the code generation
stuff?
> it
> has many convenient functions (such as the printing functions to generate
> c-code of your expressions) and one can use all other python functions. --
> and sympy is free ;)
>
> Regards,
> Vinzenz
Aaron Meurer
gundamlh
I use Sympy to compute the solution of an inverse problem (Ax = b) symbolically, where A is a 9-by-9 symbolic matrix. The computation is quite fast, however the out-printing of the solution ''x'' vector is very slow, even the first entry x[0] takes longer than 2 minutes. And the sympy.ratsimp(x[0]) is also very slow.. I am confused...
Hope for your kind response. Thanks in advance!
Hong
在 2012年3月4日星期日UTC+1下午9时49分06秒，Aaron Meurer写道：
Jason Moore
--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/65cd2d04-84ff-4eac-8a69-93c144061d83%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Richard Fateman
Aaron Meurer
printers like pprint() or sstr(), which can speed up the printing of
very large expressions.
Aaron Meurer